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waveguide phases reflecting P- and S-energy propagation. Extensive tests on envelope analysis of 
local records from different areas found that arrival times of the maxima of Pg and Lg envelopes 
increase very consistently with distance even in different tectonic regimes. Typical velocities 
being 6.1 and 3.5 km/sec respectively. These arrival time parameters are easy to extract in an 
semi-automatic manner and are highly suitable for local epicenter determinations. Extensive tests 
on locating mining explosions were conducted and on average the 'envelope* location errors j 

relative to 'true* locations were similar to those in bulletins which are based on conventional 
phase pickings. Occasionally the Pg/Pn envelope may be very weak but can be replaced by the easily | 
pickable (non-envelope) Pn-phase. Additional advantages with envelope locations are 
transportability (not overly sensitive to details of crustal structure), and that envelope 
amplitudes can be directly converted to ground motion and magnitudes. For modern stations 
envelopes can be formed in situ with low sampling rates of 1-2 Hz thus greatly reducing 
transmission costs. 

(iii) Novaya Zemlya (NZ) seismic events always stir interest among seismologists and politicians 

alike; the issue is whether a clandestine nuclear explosion has taken place at or near the former j 
USSR nuclear test site or simply an earthquake occured. In this regard the HZ event of 13 Jan, j 

1996 was of particular interest because epicenter locations as provided by various seismologist | 
groupings varied by more than 100 km and even location confidence ellipses did not overlap. We 
were asked by German colleagues to look into this epicenter location controversy and in particular 
to analyze 3-component recordings from the Norilsk station in Siberia (operated by IRIS). Other 
seismologists had ignored these recordings presumably due to * spiky* records (local electric 
outgauges) and besides relatively high noise levels. Using our new envelope processing technique 

we were able to extract Pn- (somewhat doubtful) and Sn-arrival times (reliable). Using these 
additional phase readings in combination with Pn- and Sn-phase readings from ARCESS and | 

Spitsbergen arrays we obtained an epicenter location at 75.38 N; 56.55 E that is on the NW Coast | 
of Novaya Zemlya. Part of the controversy bearing on this paxticular event was tied to an early 
epicenter location westward in the Barents Sea. We demonstrated that significant seismic phase 
information could be extracted from *poor* records and that the commonly used lASPEI travel time ■ 

tables are inadequat for Novaya Zemlya and adjacent regions. | 

(iv) An intriguing problem in seismology is the Lg-blockage across geomorphical features like > 

mountain ranges and grabens. The North Sea grabens are the classical example here with nearly \ 

total blockage for certain wave paths. Such observations are explained in terms of distortion or j 
destruction of the crustal waveguide due to crustal thinning. This hypothesis was extensively 
tested by successively thinning North Sea crystalline crustal models down to 5 km - still more 
than 50 */, of Lg-energy * survived* in the synthetic seismograms after graben passages. Very similar 
results have recently been published by other scientists albeit denoted as graben Lg-blockage. 

Even introducing as much as 8 */. RMS velocity pertubations within the graben area proper did not 
prevent Lg passage. Rough calculations indicated that a low Q-value of the order of 50 would j 

suffice for near total Lg- blockage in combination with crustal thinning effects. In interesting 
observational feature was that in cases of Lg-blockage correspondingly the Sn and its coda were 
strong and vice versa. This implies that Lg energy is leaking out of the crustal waveguide into 
the upper mantle and outside the graben area reemerge at the free surface as Sn-waves. 

(v) We are continuing our efforts on simulating seismic wave propagation in the crust and upper 
mantle. In this regard we have considered the problem of proper boundary conditions in 2D FD 
wavefield modeling. Two methods of grid boundary absorption were investigated and associated 
performances compared. Two methods of grid wave absorption are investigated and compared; the 
recent Optimal Absorbing Boundary Condition (Peng and Toksoz, J. Acoust. Soc. Am., 1994) and the 
Exponential Damping (a variant of Cerjan et. al.. Geophysics, 1985). A staggered 8th order 
accurate spatial finite-difference discretization applied to the velocity-stress formulation of 
the elastic wave equations is used. A free, plane surface is incorporated at the top boundary, 
while absorbing boundary conditions are employed along the remaining grid boundaries. We find 
conditions for achieving stable results when applying these methods in the specific computational 
environment. The respective merits are also compared using both a homogeneous model and a more 
realistic model consisting of several tilted layers. The main result is that the ED formulation 

(continued on next page) 


Unclassif _ 

SGCURUY CLASSIFICATION OF THIS PAGE 





UNCI^SSIFIHD 


ShCURlTY ClASSiFICATIONOF IHlS PAGt 


gives better edge absorbing properties than that of the OABC method. 

(vi) As mentioned in (iv) the 2D FD synthetics illustrate in an instructive manner seismic wave 

propagation in complex media albeit coda excitation appears somewhat weaJker than seen in real 
records. Intuitively, we expect that 3D models would generate relative more coda waves and in this 
regards experiments have been carried out to model scattering from free surface topography near 
the NORESS array. Three-dimensional (3“D) finite-difference (F-D) modeling of scattering from free 
surface topography has been pursued. A velocity stress formulation of the full elastic wave 
equations with exact boundary conditions for a free surface topography has been numerically 
modelled by an 8th order finite—difference method on a staggered grid. We have simulated 
scattering from teleseismic P-waves using a plane, vertically incident P-wave and real topography ; 
from a 40x40 km area centered at the NORESS array in south-eastern Norway. Snapshots and synthetic! 
seismograms of the wavefield show clear conversion from P to Rg (short period fundamental mode j 
Rayleigh) waves in an area of rough topography approximately 10 km east of NORESS. This result is j 
consistent with numerous observations. By the recent parallellization of the program code using j 
MPI (Message Passing Interface), new possibilities have been opened for modeling higher | 

frequencies and/or larger areas than so far conducted. A more complete 3-D modeling of the j 

observational axea can be realized. 

(vii) We have pursued and compared two- and three-dimensional (3-D) finite-difference (F-D) 
modeling of scattering from free surface topography. A velocity-stress formulation of the full 
elastic wave equations are combined with exact boundary conditions for the surface topography and 
numerically discretized by an 8th order F-D method on a staggered grid. We have simulated 
scattering in 2-D and 3-D from teleseismic P-waves using a plane, vertically incident P-wave and 
real topography from a 60 by 60 km area including the NORESS array in south-eastern Norway. Many 
field observations that are not easily explained by simpler 2-D cases are shown to better match 
qualitative effects from 3-D surface topography modeling. These include strong amplifications at 
hills, complex wave pattern caused by scattering, and directivity of scattered waves. Snapshots 
and seismograms show clear conversion from P- to Rg- (short period fundamental mode Rayleigh) 
waves in an area of rough topography in the vicinity of the array site. All results are consistent 
with numerous observations. By parallellization of the original software, possibilities have been 
opened for modeling with higher resolution and/or larger areas than before. 

'(viii) An essential element in nuclear test ban monitoring is that of seismic source 
classification; earthquake or underground nuclear explosion ? An additional complexity of this 
issue is that numerous signals stemming from mining and quarry explosions (chemical) are recorded 
daily and hence contribute significantly to the workload on the monitoring system. An outstanding 
feature of this kind of industrial activities is the spatial stationarity of the explosion sites 
which often are less than 1 km in aperture. In other words, a scheme for automatically recognizing 
the seismic 'signature* of chemical explosion sites would greatly reduce the daily monitoring 
workload. Basic assumption here, well established observationally, is that of seismic waveform 
similarities for closely spaced explosions and also even aftershocks (earthquakes). In many | 

industrialized countries of low seismicity more than 90 */, of seismic event recordings stem from ; 

chemical explosions and thus contribute significantly to the daily analyst workload. In this study j 
we explore the possibility of using waveforms from a priori known explosion sites (learning) for 1 
recognizing subsequent explosions from the same site excluding any analyst interference. To ensures 
high signal correlation while retaining good SNRs we used envelope transformed waveforms including- 
both the P and Lg arrivals. To ensure good spatial resolution we used multistation (network) ! 

recordings. An interpolation and approximation neural network scheme was used for learning the ’ 

computer to recognize new explosion recordings from a spesific site using detector output event 
files of waveforms only. The analysis output is a single number between 0-1 and on this scale an 
acceptance threshold of 0.4 proved appropriate. We obtained 100 */, correct decisions between two 
sets of ^site explosions' and hundreds of ^non-site' explosions/earthquakes using data files from 
the Norwegian Seismograph Network. 

(ix) Another possible approach to signal site recognition (see viii) is to replace 
multistation (z-component only) with single station 3-component recordings. Instead of using 
envelope transformed records we may introduce complex trace modulation so the signal covariance 

(continued on next page) 
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matrix would have 9 different time elements. As is well known from many source studies both high | 
and low signal frequency information contributes significantly to good classification performances 
so we formed covariance matrixes for 12 different frequency bands. In other words, a single 
3-component event recording was in our new neural net signal recognition scheme by 9 x 12 = 108 
pieces of time observation. Ground truth information was obtained for two underwater construction 
sites at Mongstad and Geiranger, W. Norway and 4-6 events were used in establishing the 
respective learning sets. Testing on hundreds of event recordings with epicenters adjacent to the 
two construction sites produced excellent results, that is all construction explosions were 
recognized as such while all other events were classified as non-site events. Probabilities of a 
’false alarm’ or a ’missed detection’ (presuming score values to be Gaussian) were of the order of 
1:1,000,000. We also checked the relative merits (contributions) of individual filter bands and found 
that the 3 - 6 Hz and the 8 - 12 Hz segments were most informative. This is somewhat surprising since 
these two signal segments are also most informative in the context of discrimination between 
explosion and earthquake sources. [ 
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1 Objectives 

• Foster a better understanding of seismic wave propagation in complex lithospheric media. 

• Principal modeling tool is 2-D and 3-D finite difference (FD) synthetics for models also 
including surface topography. 

• Corroborating synthetics through analysis of array and network data. 

• Use ’synthetic’ knowledge in design of event classification schemes. 

• Seismic network performance and real time event location. 

• Automated explosion site spesific signal recognition 


2 Status of effort 

Much efforts have been spent on expanding to 3D models and streamlining our finite difference 
(FD) wavefield simulation scheme including code adaption for parallel computers. The latest 
extension here is to include intrisic attenuation or damping that is numeric solutions of the vico- 
elastic wave equations in 2D and 3D. The 3D FD modeling of topography wavefield responses in the 
NORESS array siting area are rated a success; we are able to reconstruct realistically the P-to-Rg 
scattering fields with a 'secondary’ source in the Bronkeberget hills. Using the 2D FD variant we 
have also simulated T^f-blockage across the North Sea graben structures. Partial blockage of Lg is 
tied to crustal thinning with Lg leaking out of the waveguide in the form of 5n-waves. Significant 
Lg decay is also caused by Lg-io-Rg conversions in the overlying sediment strata and by intrinsic 
attenuation. In a CTBT monitoring context, near real time event classification is essential at 
least in the local epicenter distance range. • First step of our efforts here have been to design an 
automatic event detection and epicenter location scheme. It is very robust and when the locations 
are compared with known mining sites, the results appear very accurate. More realistic real time 
testing is planned in cooperation with colleagues at the German National Data Center (NDC) in 
Hannover. Jointly with our German colleagues we analyzed a presumed ’suspecious’ event taking 
place on the NW Coast of Novaya Zemlya 13 Jan. 1996; outcome was that there was no strong 
evidence that this event could be a clandestine nuclear explosion. Numerous mining and quarry 
sites generate seismic signals which constitute a heavy workload on any CTBT monitoring system. 
Since explosion sites are spatially stationary signal records from a spesific site are rather similar. 
We have exploited this this observational fact by designing a neural net scheme for automatically 
recognizing site spesific signals thereby eliminating any need for analyst interference. In fact, 
we have exploited two approaches here; one using signal envelopes from a set (3 - 6) of network 
stations (fiexibel choices) and the other one tied to single station 3-component recordings. In the 
latter case we used the ’complex modulated’ covariance matrix for a 12-suit of bandpass filters so 
one event record was equivalent to 9x12 = 108 pieces of time information. Use of single station 
recordings were in part motivated by network timing errors. Both ’site recognition’ schemes had 
excellent performances when tested on local explosion recordings in W. Norway for which ’Ground 
Truth’ information was available. 
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Final work effort - a personal farewell 


The Principal Investigator, Eystein S. Husebye, was introduced to US Nuclear Monitoring Research 
(VELA Uniform Program) a long time ago, that is, in May 1966 when offered a Post. Doc. 
Fellowship for studies at MIT, Cambridge, MA. Afterwards employed as Director of Research 
(1968 - 1993) at NORSAR (The Norwegian Seismic Array), Kjeller, Norway which construction 
and subsequent operation for many, many years were mainly sponsored by AFTAC and ARPA. 
From 1991 - 1997 PFs (our) seismic monitoring research have been sponsored by the Air Force 
Office of Scientific Research (AFOSR) Bolling AFB and managed by Phillips Laboratory, Hanscom 
AFB. Pi’s longstanding involvement with various US Defense and Air Force agencies have been 
intellectually very stimulating and personally rewarding in terms of an above average scientific 
and seismological career. In the latter case, to be more spesific, editor of 2 books on nuclear 
test bam monitoring, more than 150 technical reports and scientific journal contributions and 
organizer/convener for numerous international symposia, conferences and NATO ASIs. In short, 
Eystein S. Husebye is very grateful for the many research opportunities offered by participation 
in the many US sponsored scientific programs within the fields of nuclear test ban monitoring and 
basic seismology these wordings are just a very modest way of expressing this gratitude. 
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4 Accomplishments 


In the following a short summary of each subject is given. The complete papers appear in appen¬ 
dices, Some of the topics listed were included in our previous Interim Project Report; reason for 
this is that the relevant reports (4.4, 4.5, 4.6) have been significantly modified or is now appear¬ 
ing as journal publication. Sec. 4.1 is an exception since organizing a NATO Advanced Study 
Institute and subsequently editing an important book on Nuclear Monitoring is rated a major 
accomplishment and hence repeated here. 

4,1 NATO ASI symposium &: book on CTBT monitoring 

Our NATO ASI Symposium in Alvor, Portugal 23 Jan-1 Feb, 1995 on CTBT monitoring was 
rated a great success. The final outcome of this undertaking was the book: E.S. Husebye and A.M. 
Dainty (eds): Monitoring a Comprehensive Test Ban Treaty, Kluwer Academic Publ., Doordrecht, 
The Netherlands, pp 836, 1995. The book contains contributions from many scholars on most 
CTBT issues albeit with emphases on seismic monitoring. More than 100 key lecturers and 
students attended the symposiom including 10 scientists and nuclear physicists from the former 
Soviet Union. 


4.2 Local event location in an automatic manner 

Local event recordings are complex in the sense that relevant P- and S-phases vary in an unpre¬ 
dictable manner even between closely spaced stations; thus, manual analysis of such records is still 
commonplace. Our approach to solving this long standing problem of observational seismology is 
to bandpass filter (3-6 Hz) to ensure good SNR and then form envelopes to ensure simple signals 
across a seismograph network. The physical basis is that Pg and Lg are crustal waveguide phases 
reflecting P- and S-energy propagation. Extensive tests on envelope analysis of local records from 
different areas found that arrival times of the maxima of Pg and Lg envelopes increase very con¬ 
sistently with distance even in different tectonic regimes. Typical velocities being 6.1 and 3.5 
km/sec respectively. These arrival time parameters are easy to extract in an semi-automatic man¬ 
ner and are highly suitable for local epicenter determinations. Extensive tests on locating mining 
explosions were conducted and on average the ‘envelope' location errors relative to ‘true' locations 
were similar to those in bulletins which are based on conventional phase pickings. Occasionally 
the Pg/Pn envelope may be very weak but can be replaced by the easily pickable (non-envelope) 
Pn-phase. Additional advantages with envelope locations are transportability (not overly sensitive 
to details of crustal structure), and that envelope amplitudes can be directly converted to ground 
motion and magnitudes. For modern stations envelopes can be formed in situ with low sampling 
rates of 1-2 Hz thus greatly reducing transmission costs. 

4.3 Relocating the 13 January 1996 Novaya Zemlya seismic event 

Novaya Zemlya (NZ) seismic events always stir interest among seismologists and politicians alike; 
the issue is whether a clandestine nuclear explosion has taken place at or near the former USSR 
nuclear test site or simply an earthquake occured. In this regard the NZ event of 13 Jan, 1996 was 
of particular interest because epicenter locations as provided by various seismologist groupings 
varied by more than 100 km and even location confidence ellipses did not overlap. We were asked 
by German colleagues to look into this epicenter location controversy and in particular to analyze 3- 
component recordings from the Norilsk station in Siberia (operated by IRIS). Other seismologists 
had ignored these recordings presumably due to 'spiky' records (local electric outgauges) and 
besides relatively high noise levels. Using our new envelope processing technique we were able 
to extract Pn- (somewhat doubtful) and Sn-arrival times (reliable). Using these additional phase 
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readings in combination with Pn- and Sn-phase readings from ARCESS and Spitsbergen arrays 
we obtained an epicenter location at 75.38 N; 56.55 E that is on the NW Coast of Novaya Zemlya, 
Part of the controversy bearing on this particular event was tied to an early epicenter location 
westward in the Barents Sea. We demonstrated that significant seismic phase information could 
be extracted from ’poor’ records and that the commonly used lASPEI travel time tables are 
inadequat for Novaya Zemlya and adjacent regions. 


4.4 X^f-blockage across North Sea grabens 

The North Sea T^^-blockage for wave paths across the crustal graben structures is a well-established 
observational fact. Analysis of such observations implies that Z-p-blockage takes place in graben 
areas associated with sedimentary basin formation and crustal thinning. These intriguing obser¬ 
vations have triggered many theoretical studies aimed at highlighting specific Lg loss mechanisms 
albeit so far with only moderate success. Our approach to this problem is to simulate seismic 
wavefield propagation through the crustal waveguide using 2D finite difference techniques. ^From 
oil exploration works in the North Sea the graben structures are known in detail which enabled us 
to use realistic crustal models in our Lg synthetics. In the most extreme model tested, the crys¬ 
talline crust thickness beneath the graben amounted to only 5 km while the overlying sedimentary 
pile is nearly 10 km thick. At the base of the crust the Moho is elevated nearly 10 km below the 
graben. This model has similarities to the oceanic crustal waveguide where total T^^-blockage is 
claimed for path lengths exceeding 100 km. The synthetic wavefields are displayed in terms of 
snapshots, semblance velocity analysis and time-space RMS amplitudes. The dominant structural 
Lg loss mechanisms are Lg-to~Rg conversions (scattering) in the sediments and S-wave leakage out 
of the crustal waveguide and into the upper mantle. Part of these upper mantle .S-waves return to 
the crust and appears as Sn coda. Observationally, strong 5n-phases of long duration are often 
associated with weak Ty-phases and vice versa. Our synthetics produced Lg amplitude decay 
amounting at most to 6-10 dB while observational data imply blockage amounting to 15-20 dB. 
The latter is equivalent to Pn-Lg magnitude difference of nearly one magnitude unit. The main 
outcome of this study is therefore that Lg-wa,ye propagation is very robust and that a dominant 
blockage effect associated with intrinsic attenuation, that is Q values of the order of 100 at 2 Hz 
for a path length of minimum 100 km, is necessary to conform to observations. 


4.5 Boundary conditions in 2D FD modelling 

Two methods of grid wave absorption are investigated and compared; the recent Optimal Ab¬ 
sorbing Boundary Condition (Peng and Toks0z, J. Acous. Soc. Am., 1994) and the Exponential 
Damping, a variant of Cerjan et. al. (Geophysics, 1985). A staggered 8th order accurate spatial 
finite-difference discretization applied to the velocity-stress formulation of the elastic wave equa¬ 
tions is used. A free, plane surface is incorporated at the top boundary, while absorbing boundary 
conditions are employed along the remaining grid boundaries. We find conditions, for achieving 
stable results when applying these methods in the specific computational environment. The re¬ 
spective merits are also compared using both a homogeneous model and a more realistic model 
consisting of several tilted layers. 

4.6 3D FD modeling of scattering from topography 

Three-dimensional (3-D) finite-difference (F-D) modeling of scattering from free surface topog¬ 
raphy has been pursued. A velocity stress formulation of the full elastic wave equations with 
exact boundary conditions for a free surface topography has been numerically modelled by an 8th 
order finite-difference method on a staggered grid. We have simulated scattering from teleseismic 
P-waves using a plane, vertically incident P-wave and real topography from a 40 x 40 km area 
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centered at the NORESS array in south-eastern Norway. Snapshots and synthetic seismograms 
of the wavefield show clear conversion from P to Rg (short period fundamental mode Rayleigh) 
waves in an area of rough topography approximately 10 km east of NORESS. This result is con¬ 
sistent with numerous observations. By the recent parallellization of the program code using MPI 
(Message Passing Interface), new possibilities have been opened for modeling higher frequencies 
and/or larger areas than so far conducted. A more complete 3-D modeling of the observational 
area can be realized. 

4.7 3-D versus 2-D finite difference seismic synthetics 

We have pursued and compared two- and three-dimensional (3-D) finite-difference (F-D) modeling 
of scattering from free surface topography. A velocity-stress formulation of the full elastic wave 
equations are combined with exact boundary conditions for the surface topography and numeri¬ 
cally discretized by an 8th order F-D method on a staggered grid. We have simulated scattering 
in 2-D and 3-D from teleseismic P-waves using a plane, vertically incident P-wave and real to¬ 
pography from a 60 by 60 km area including the NORESS array in south-eastern Norway. Many 
field observations that are not easily explained by simpler 2-D cases are shown to better match 
qualitative effects from 3-D surface topography modeling. These include strong amplifications at 
hills, complex wave pattern caused by scattering, and directivity of scattered waves. Snapshots 
and seismograms show clear conversion from P- to Rg- (short period fundamental mode Rayleigh) 
waves in an area of rough topography in the vicinity of the array site. All results are consistent 
with numerous observations. By parallellization of the original software, possibilities have been 
opened for modeling with higher resolution and/or larger areas than before. 

4.8 Recognizing explosion sites without seismogram readings 

An essential element in nuclear test ban monitoring is that of seismic source classification; earth¬ 
quake or underground nuclear explosion ? An additional complexity of this issue is that numerous 
signals stemming from mining and quarry explosions (chemical) are recorded daily and hence 
contribute significantly to the workload on the monitoring system. An outstanding feature of 
this kind of industrial activities is the spatial stationarity of the explosion sites which often are 
less than 1 km in aperture. In other words, a scheme for automatically recognizing the seismic 
’signature' of chemical explosion sites would greatly reduce the daily monitoring workload. Ba¬ 
sic assumption here, well established observation ally, is that of seismic waveform similarities for 
closely spaced explosions and also even aftershocks (earthquakes). In many industrialized coun¬ 
tries of low seismicity more than 90 % of seismic event recordings stem from chemical explosions 
and thus contribute significantly to the daily analyst workload. In this study we explore the possi¬ 
bility of using waveforms from a priori known explosion sites (learning) for recognizing subsequent 
explosions from the same site excluding any analyst interference. To ensure high signal correlation 
while retaining good SNRs we used envelope transformed waveforms including both the P and 
Lg arrivals. To ensure good spatial resolution we used multistation (network) recordings. An 
interpolation and approximation neural network scheme was used for learning the computer to 
recognize new explosion recordings from a spesific site using detector output event files of wave¬ 
forms only. The analysis output is a single number between 0-1 and on this scale an acceptance 
threshold of 0.4 proved appropriate. We obtained 100 % correct decisions between two sets of ‘site 
explosions’ and hundreds of ‘non-site’ explosions/earthquakes u^^^ files from the Norwegian 

Seismograph Network. 


6 



4.9 Signal site recognition using single station 3-component records 

As mentioned, seismic waveform similarities for closely spaced explosions are well established ob¬ 
servation ally. Also, encouraged by our previous site recognition results we considered additional 
flexibility in our approach to this problem essentially how to incorporate additional seismic record 
information in the neural network analysis. Another motivation here was to safeguard against 
technical defects like timing errors which otherwise would ruin multistation analysis. The alter¬ 
native strategy choosen for extracting relative comprehensive signal attributes was that of using 
single station 3-component (3C) records. By introducing complex demodulated record operations 
(including both the P- and S-wave arrivals and part of the coda) the 3C covariance matrix would 
contain 3x3 elements of different time information. From other source classification studies we 
know that both explosion and earthquake spesific information are contained in several frequency 
bands so we formed covariance matrixes for a suit of 12 different frequency bands. In other words, 
a single 3C station event record produced 9x12 = 108 pieces of time information. Having access 
to ’ground truth’ information from 2 underwater construction sites (Mongstad and Geiranger) 
we had 4-6 events (enough) for site recognition learning sets. Obtained results were excellent; 
all events from the construction sites (not in learning sets) were classified as such while all other 
events (more than 200) were classified as non-site events whether source was an explosion or an 
earthquake. Presuming log-normal distribution of the neural network output score factor we esti¬ 
mated that the probabilities of false alarm or missed ’detection’ were less than 1:1,000,000. Also, 
spatial resolution was good (less than 10 km) even perpendicularly to the station-site azimuth 
direction. During this 3C site recognition exercise we found in the epicenter map used, 2 cluster of 
events with which origin times were mainly between 10. - 14. hours that is prime working hours. 
Selecting events directly from bulletin for learning sets but without ground truth information at 
hand did produce results somewhat less convincing in comparision to those obtained for Mongstad 
and Geiranger. This kind of problems, lack of ground truth information for appearant event clus¬ 
ters, would in future be addressed through an initial multivariate analysis scheme in order to sort 
out ’most similar’ event recordings and hence the likely event candidates for learning sets. 
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PREFACE 


This NATO Advanced Study Institute addressed the problem of verifying compliance by 
signatories of a Comprehensive Test Ban Treaty (CTBT), including setting up a 
monitoring system that will detect treaty violations. The verification problem 
comprises a blend of political, scientific and technical issues, which have often been 
treated separately. The requirements of the system are set by political considerations; it 
must be based on scientifically sound principles; and it must be implemented 
technically in the real world. 

The issue of a Comprehensive Test Ban Treaty (CTBT) has been on the U.N. 
political agenda for almost 40 years. As this is being written, it is an active topic: the 
Conference on Disarmament (CD) is working to draft such a treaty which will then be 
made available for any nation that wishes to sign. Another aspect of a CTBT is its 
place in the larger issue of non-proliferation of nuclear weapons. Again, this is an 
active topic as this is being written: negotiation of a CTBT was explicitly stated as a 
goal when the Non-Proliferation Treaty was recently renewed indefinitely. In the 
technical field, the Group of Scientific Experts (GSE), convened by the CD, is 
conducting its third Technical Test, known as GSETT-3. The intention of GSETT-3 is 
to produce a prototype system for seismological monitoring that, among othCT goals, 
can serve as a practical example for the negotiating community of what is technically 
feasible. Also, the CD has established special working groups on the non-seismic 
techniques of hydroacoustic, radionuclide, and mfiasound monitoring, and on on-site 
inspection methods, as part of the process of drafting a CTBT. 

In these Proceedings, the main focus will be on the scientific and technical aspects 
of seismic monitoring of underground nuclear tests. This particular topic was chosen 
because it is one of the main technical challenges o^monitoiing. However, we have 
also included the political issues and non-seismiC- techniques. Accordingly, the 
Proceedings of the ASI have been arranged so that a general overview of the problem of 
monitoring is given in the first section, with a more detailed exposition of the technical 
practice and problems of seismic monitoring foUowing, and papers on future directions 
in the last section. The overview section includes contributions about political 
background and developments, and technical reviews of the latest practice in testing and 
the non-seisnuc monitoring methods, including on-site inspection. Some, but not all, 
of the non-seismic methods are used for explosions which are not contained. These 
must be monitored, but it must be bom in mind that current testing practice is to 
contain tests underground, and any compliance system must contain provisions fcr 
monitoring such events; hence, the focus on seismic monitoring in this ASI. 

The specific nature of a CTBT, as compared to other treaties such as the Limited 
Test Ban Treaty banning all tests except underground or the Threshold Test Ban Treaty 
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' banning underground tests above a certain rather large threshold, imparts its own unique 
characteristics to a seismic monitoring system. First, a CTBT is intended to be 
multilateral and worldwide rather than bilateral or trilateral, thus any monitoring system 
must cover the globe and all environments: underground, underwater, atmosphere and 
space. Second, and technically much more important, all tests of any size would be 
banned under a CTBT, including small events. To monitor small underground events, 
it is necessary to record seismrigrams at relatively short ranges of 100-2000 km from 
the event, known in seismology as the regional distance range. This is a major change 
from current operational procedure in monitoring which is based at monitoring at 
distances beyond 2000 km, the teleseismic distance range: it has been known since the 
earliest days of seismology that seismograms at regional distances have quite different 
characteristics and frequency content than seismograms at teleseismic distances, due to 
the different parts of the earth that are sampled during propagation of the signal between 
source and receiver. 

The new technologic^ developments required by this shift to regional monitoring 
- form an important part of the papers presented here. The operational process of seismic 
monitoring can be divided into detection of all seismic signals, location of all events, 
identification of the seismic events which are underground nuclear tests (known as 
’^discrimination” in seismic monitoring), and estimation of the size (yield) of the 
nuclear explosions. In assessing the new problems of regional monitoring, it must of 
course be remembered that expertise from other fields in seismology, especially 
earthquake recording, is available, and indeed many of the resources of knowledge and 
equipment already in place will surely be used in monitoring a CTBT. However, such 
knowledge and equipment has been acquired to meet the needs of these fields, and will 
not always fit the requirements of seismic monitoring. The most pressing case is 
discrimination, since this is a subject that is umque to seismic monitoring. 
Furthermore, because of the higher frequency of regional signals compared to 
teleseismic, the most effective means of discrimination discovered for teleseismic 
signals (MS:mb) cannot be used for small events and new methods must be found. 
Accordingly, particular attention has been paid to discrimination. 

The structure of the Proceedings is designed to cqver the problems listed above. 
“After opening with an overview section on Principal Political and Techmeal Test Ban 
Issues, a section on Monitoring Technologies gives a summary of all the methods 
being considered for monitoring a CTBT. Following this, the section Explosion and 
Earthquake Source Modelling examines underground nuclear tests and other disturbances 
as seismic sources. This examination includes the problem of decoupling (i.e., 
decoupling the explosion from the solid earth by detonating it in a cavity), whiefa^ is 
probably the most credible way in which monitoring might be defeated. Following 
this the problem of monitoring such seismic sources is discussed in some detail, with 
sections on Seismic Networks, including detection and location of seismic events, 
Signal.Processing and Seismic Wave fiopagation, coy speci fic 

to the seismic arrays used in regional monitoring and propagation of seismic waves at 
regional distances; and Seismic Source Discrimination Technologies, especially at 
regional distances. These topics cover three of the four classical steps of monitoring. 
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namely, detection, location, and discrimination; the signal analysis and wave 
propagation are included as the scientific basis for accomplishing this. A review -papa 
on existing methods of yield estimation is included in the section on source modelling, 
but this topic has not been treated at greater length because in a CTBT, all tests of any 
size win be proscribed. 

The subject of monitoring underground nuclear explosions has its points of 
controversy, technical and otherwise. In the NATO ASI, and in these Proceedings, our 
position has been to present as many points of view as possible and to encourage the 
contributors to put down there own points of view. We have made no effort to 
reconcile differing points of view, but simply present them for the reader to make his or 
her own judgments. Within these constraints, we hope these Proceedings give an 
accurate and succinct review of the state-of-the-art (or state-of-the-argument) of 
Monitoring a Comprehensive Test Ban Treaty. It can be said that the practicaUy 
unanimous opinion of the ASI participants, as expressed in the scheduled discussion 
sessions, is that a Comprehensive Test Ban Treaty can be monitored in the technical 
sense, at least down to some yield; the arguments are over what that yield is. 

A great deal of effort goes into organizing and running an ASI. We would like to 
express our sincere thanks to all who worked, often late into the night for man y nights, 
to make this Conference possible. Particular thanks go to Prof. F. Tomas, the State 
Secretary of Science and Technology of Portugal, and Prof. R. Ribiero, President of 
Junta Nacional de Investigacao Cientifica e Tecnologica, Lisboa, Portugal, who found 
time in their busy schedules to attend and speak at the Opening Ceremonies. We are 
also indebted to our fellow Organizing Committee members Prof- Alexander Kushnir of 
the International Institute of Earthquake Prediction Theory aixl Mathematical 
Geophysics, Moscow, Russia, and Prof. Luis Mendes-Victor of the University of 
Lisboa, Portugal. Nothing would have been possible without the work of the Local 
Organuing Committee consisting of Prof. Mendes-Victor, Prof. Paula Costa, Prof. 
Carlos OUveira, Dr. Luisa Senos, Dr. Luis Matia, and Fernanda Dias, aU of the 
University of Lisboa; they have our gratitude for a smoothly run and enjoyable 
meeting. 

-tfr 

August 1995 


Eystem Husebye 
Bergen, Norway 


Anton Dainty 
Boston, U.S.A.. 
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Robust and Reliable Epicenter Determinations: 
Envelope Processing of Local Network Data 

by Eystein S. Husebye, Bent O. Ruud, and Anton M. Dainty* 

Abstract Local event recordings are complex in the sense that relevant P and S 
phases vary in an unpredictable manner even between closely spaced stations; thus, 
manual analysis of such records is still commonplace. Our approach to solving this 
long-standing problem of observational seismology is to bandpass filter (3 to 6 Hz) 
to ensure good signal-to-noise ratio SNR and then form envelopes to ensure simple 
signals across a seismograph network. The physical basis is that Pg and Lg are crustal 
wave-guide phases reflecting P- and 5-energy propagation. Extensive tests on en¬ 
velope analysis of local records from different areas found that amval times of the 
maxima of Pg and Lg envelopes increase very consistently with distance even in 
different tectonic regimes, typical velocities being 6,1 and 3.5 km/sec, respectively. 
These arrival-time parameters are easy to extract in a semi-automatic manner and 
^ are highly suitable for local epicenter determinations. Extensive tests on locating 

mining explosions were conducted, and on average, the “envelope” location errors 
relative to “true” locations were similar to those in bulletins that are based on 
conventional phase pickings. Occasionally, the PgIPn envelope may be very weak 
but can be replaced by the easily pickable (nonenvelope) Pn phase. Additional ad¬ 
vantages with envelope locations are transportability (not overly sensitive to details 
of crustal structure), and that envelope amplitudes can be directly converted to ground 
motion and magnitudes. For modem stations, envelopes can be formed in situ with 
low sampling rates of 1 to 2 Hz, thus greatly reducing transmission costs. 


Introduction 

A common assumption in local seismogram analysis is 
that the records are dominated by the four phases Pn, Pg, 
Sn, and Lg (replaced by Sg at short distances). These phases 
are associated with specific travel-time curves and propa¬ 
gation paths that in turn are used in epicenter determinations, 
ftp This simple approach is now being questioned. For example, 
Pg and Lg are not adequately described as simple rays; they 
are essentially wave-guide phases including scattering con¬ 
tributions from crustal heterogeneity (Hestholm etaL, 1994.) 
At distances beyond about 100 km, the direct Sg signal is 
rather weak, hence the change of notation to Lg —the dom¬ 
inant crustal wave observable occasionally out to 3000 km 
(Nuttli, 1973; Hansen et aL, 1990). Since normally Pg and 
Lg are dominant over Pn and Sn, the signal envelopes of P- 
and 5-wave arrivals are indicative of P- and 5-energy trans¬ 
port in the crustal wave guide. An indication for this is that 
the L^-wave group velocity is reported in the interval 3.4 to 
3.6 km/sec for large parts of the continents (Mykkeltveit and 
Husebye, 1981). 
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For local epicenter determinations, preferences are to 
work with high signal frequencies {f>2 Hz) where signal- 
to-noise ratio (SNR) is relatively high, thereby accepting the 
penalty of complex signal records. Low-frequency filtering 
typically results in simple records but poor SNR, thus this 
option is not of much interest. The greatest problem is that 
the shear wave rain {Sn and/or Lg) that typically is the most 
prominent feature in the seismogram often has an emergent 
onset. Hence, reliable phase picks are difficult to accomplish 
automatically, so these efforts remain essentially manual. 
When we contemplated this problem, the envelope process¬ 
ing introduced in NORSAR array surveillance appeared at¬ 
tractive (Ringdal et al, 1975; Nuttli, 1973): Could this be a 
successful approach for automated analysis of complex local 
records as well? By envelope transforming such records, we 
obtain traces with excellent SNR (filter band 3 to 6 Hz used), 
but signals are not identical since local network apertures 
are often large. However, both PG and LG (envelope nota¬ 
tions for Pg and Lg) peaks are sufficiently similarly shaped 
to warrant combined analysis operations based on ID signal 
attributes; that is, envelope parameter variations depend on 
epicenter distance alone; the original signal waveforms do 
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not. By plotting envelope traces (Fig. 2), it is easy to verify 
the validity of this assumption (Fig. 3), and besides, PG and 
LG (their envelope peaks) propagate with velocities of ap¬ 
proximately 6.1 and 3.5 km/sec, respectively, as expected 
for crustal wave-guide propagation. Ryzhikov et ai (1996) 
designed an ingenious scheme for joint PG and LG power 
location in time and geographical space, equivalent to epi¬ 
center determination. Conceptually, this has some similari¬ 
ties to global beamforming and threshold monitoring (Ser- 
eno and Bratt, 1989; Kvsema and Ringdal, 1995). In essence, 
envelope location may be considered a variant of/-/: anal¬ 
ysis, tomography backprojection, or a “brute force” search 
for maximum power in geographical space as a function of 
time. 

There are drawbacks with “envelope power” location 
schemes, namely, the common necessity of weighting the 
dominant LG signal, lack of resolution due to a broad LG 
peak, and loss of the relatively weak PG beyond about a 
400-km distance. There is also the well-known problem that 
for P-only observations for events outside a network, gross 
errors may occur because origin time and distance are not 
well resolved, which is an incentive to use both P and S 
arrivals in spite of the difficulties. To avoid weighting and 
related problems, we decided to use the easily pickable PG 
and LG peak arrival times in a conventional epicenter lo¬ 
cation scheme as demonstrated below. 

Time-Domain Envelope Location 

Most local networks operate in a trigger mode; that is, 
only segmented data containing real signals are transmitted 
to the local hub for further analysis. Such data are used here, 
and after prehltering in the 2 to 4 or 3 to 6 Hz passband, 
envelopes are formed: LG peaks are picked easily by com¬ 
puter as the maximum in the segment. A preliminary loca¬ 
tion is then made from these LG readings that in turn is used 
to predict the PG arrival windows, with subsequent PG ar¬ 
rival-time picking. In cases where PG readings are problem¬ 
atic, we may instead aim at reading the first phase arrival, 
Pg or Pn (Ruud and Husebye, 1992). Then with PG and LG 
(or Pn and LG, Pg and LG) paired arrival times, reliable 
epicenter determinations are feasible using any of several 
location algorithms. Our preference is the Ruud et al (1993) 
grid-search routine, but we also used the fast Interval Arith¬ 
metics scheme of Tarvainen et al. (1997). Envelope traces 
for an event detected by the local Norwegian Seismic Net¬ 
work (NSN, Fig. 1) are shown in Figure 2, illustrating the 
ease by which PG and LG can be picked. 

In order to find appropriate travel-time curves for the 
PG and LG peaks, we proceeded in the following manner. 
First Pg arrival times (first onsets) for nearby stations (dis¬ 
tance-less than 100 km> were read.'Since- station and epi¬ 
center coordinates are known, accurate origin times were 
easy to obtain assuming a Pg velocity of 6.1 km/sec. Then 
by plotting travel times for PG and LG peaks, it was simple 
to calculate linear travel-time curves in a least-squares sence 
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0‘ 10* 20* 30* 



0* 10* 20* 30* 


Figure 1. The Norwegian Seismograph Network 
and Scandinavian array stations used in analysis. Ex¬ 
plosion sites for which envelope location accuracy 
were tested (TI = Titania, GE = Geiranger, RA = 
Rana, SV = Svartisen) and an earthquake (EQ) lo¬ 
cation examined using the arrays are also shown. The 
Titania and Geiranger explosions are only recorded 
by S. Norway stations while the Rana and Svartisen 
explosions are only recorded by stations north of 
65" N. 


(Fig. 3) for Titania and Geiranger explosion recordings. The 
scattering in PG arrival times is larger than that of LG —not 
unexpected since the crustal wave guide is more pronounced 
for shear waves. Further refinements of these ID time curves 
do not seem practical in view of lateral velocity and thick¬ 
ness variations of the crust (e.g., see Kinck et al, 1993). A 
preferable alternative is to introduce spatial time corrections 
here (Shearer, 1997). 

Network Analysis and Location Results 

Mining Explosions in Norway 

We used local explosion records from NSN stations 
(Fig. 1) to test the robustness and accuracy of our novel 
envelope epicenter determination scheme. Location perfor¬ 
mances are shown in Figure 4 for the four known explosion 
sites Titania, Geiranger, Rana, and Svartisen (Fig. 1). For 
Geiranger locations, station MOL is rather critical at an epi¬ 
center distance of only 55 km. However, at such close 
ranges, PG and LG are not well separated, so for MOL, PG 
time picks are available for only about 30% of the events. 
Geiranger locations by which all PG time picks (including 
the missing ones for MOL) are replaced by those of the an- 
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Figure 2. (a) Record displays for an explosion in the Titania quarry in southern 

Norway recorded by the Norwegian Seismograph Network stations. A 3 to 6 Hz band¬ 
pass filter is used to enhance SNR. (b) Envelope transforms after smoothing of the 
records in (a), including group velocity markers (in km/sec) tied to the focal parameters 
displayed at the top of the figure. Epicenter azimuth and distance are listed below the 
station codes. 


alyst (that is Pg and Pn arrival times) are also shown in 
Figure 4. Our envelope locations are occasionally even 
closer to the true ones than those appearing in the local bul¬ 
letins. The envelope epicenters are based on NSN station 
life records only while the bulletin solutions incorporate data 
from arrays to the east. This should give relatively better 
east-west resolution in the bulletin solutions, but this is not 
fully realized in practice due to inadequate travel-time tables. 
In Table 1, we give statistics of the location errors for both 
bulletin and envelope epicenter estimates. 

Fennoscandian Array Records 

Records from array center sensors of the arrays that 
have been offered for the International Monitoring System 
(Anonymous, 1996) are displayed in Figure 5 for an earth¬ 
quake north of Bothnian Bay (Fig. 1). Again, predicted FG 
and LG arrival times coincide well with respective record 
peaks, although PG is weak beyond 400 km. As noted 
above, for location purposes, we may replace PG with Pn. 
Since we pick arrival times tied to wavelet peaks, the SNR 


differences between an array beam trace and a single sensor 
envelope is moderate. 

Envelope Analysis of Other Network Recordings 

Since Pg, Pn, Lg, and Sn are globally observed crustal 
phases, we intuitively expect that envelope processing and 
subsequent epicenter determinations as described earlier are 
transportable to any seismic network area. The validity of 
this hypothesis has been tested on network recordings from 
Israel (Shapira et al, 1996; Shoubik et al, 1996), Italy and 
Germany (Ryzihkov et al., 1996), and also New England 
(ourselves). The outcome of all these studies were similar to 
those obtained here, namely, that envelope records are dom¬ 
inated by PG and in particular LG phases having linear 
travel-time curves with velocities of approximately 6.1 and 
3.5 km/sec, respectively. In the extreme, original waveform 
-records were too complex for reading more than onset time 
of first-arriving P phase, while after envelope transforma¬ 
tion, both PG and LG were easily pickable. In other words, 
envelope processing is clearly a viable alternative to con¬ 
ventional analysis of waveform recordings. 
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Figure 3. PG and LG peak amplitude travel times 
as a function of distance. These curves were obtained 
from observed PG and LG arrival times for Titania 
and Geiranger explosions. Since site coordinates are 
exactly known, accurate origin times were calculated 
from Pg arrival times (for distances less than 100 km) 
with a presumed Pg-phase velocity of 6.1 km/sec. 

. There are some scatters in the observations (system 
timing errors?), but the fitted straight lines have rea¬ 
sonable slopes. For even better fit to given travel-time 
curves, station corrections should be introduced. 


Discussion and Conclusion 

In this study, we have demonstrated that simple enve¬ 
lope processing is well suited for analysis of complex, local 
network recordings. The outstanding feature is the ease by 
which a substitute for the Lg phase can be picked. This phase 
is problematic in conventional analysis schemes whether 
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working with single-component, three-component, or array 
records. The envelope records give excellent SNRs; picking 
arrival times at maximum wavelet amplitude means that 
meaningful signal attributes can be extracted from single¬ 
station records that otherwise may require array recordings. 
Envelope records from a seismograph network can be used 
for real-time monitoring purposes in much the same way as 
arrays (Ryzhikov et al, 1996; Shapira a/., 1996; Shoubik 
et aL, 1996). Magnitude determinations are easily performed 
from maximum amplitudes of PG and LG wavelets, as dem¬ 
onstrated by Mendi and Husebye (1994,1995). Furthermore, 
envelope transform operations can easily be performed in 
the field, thus reducing transmission loads significantly in 
view of envelope sampling rates of 1 to 2 Hz compared to 
40 to 100 Hz for waveforms. 

Critical for accurate locations of local events is use of 
proper travel-time curves. In Fennoscandia, crustal thick¬ 
nesses vary from 20 km (offshore Norway) to about 50 km 
in central Finland (Kinck et aL, 1993). All the seimologicai 
centers in the region use travel-time curves calculated from 
two- and three-layer ID models with crustal thicknesses 
varying from 30 to 45 km. This may explain why gross 
location errors occasionally occur—in particular, in combi¬ 
nation with poor network configuration like for events oc¬ 
curring outside the network. As discussed earlier, travel-time 
curves for PG and LG are much more easily established in 
comparison to those of the conventional crustal phases. To 
our surprise, occasionally, some gross mislocations were due 
to timing errors at one or more NSN stations despite on-line 
access to the GPS (satelite) clock. The reason for this appears 
to be the software used for calibrating the internal clock 
versus the GPS clock—apparently not an uncommon prob¬ 
lem in local network operations. 

What are the potential drawbacks with envelope pro¬ 
cessing? It is not foolproof; PG may be larger than LG (es¬ 
pecially at distances less than 100 km), while Lg blockage 
is common for offshore events (Mendi et aL, 1997). These 
kinds of exceptional wave-propagation features would also 
cause problems in conventional analysis, emphasizing the 
need for analyst screening of automated bulletin results. Fo¬ 
cal depths cannot be estimated from envelope wavelet at¬ 
tributes since these parameters are not particularly sensitive 
to source depth. Perhaps the most severe objection to en¬ 
velope processing came from a colleague claiming that it 
“takes the fun out of seismology.*’ 

Envelopes of different events stemming from the same 
explosion site appear to be identical in form, implying that 
for such recordings there should be no need to undertake any 
form of phase reading. We are now testing a neural network 
scheme aimed at teaching the computer to recognize site- 
specific explosion recordings (Fedorenko et (3/.,^f99X)r-thus_ 
attempting to mimic seismogram interpretation skills of a 
trained analyst (Israelsson, 1990; Riviere-Barbier and Grant, 
1993). 


pg4 


























/ssa/88U97007 1 2/16/97 04:28PM Plate # 0 


pg5 


Short Notes 





(a) 



(c) 




Figure 4. (a) Envelope location experiment using Titania mining explosions. 

Southern Norway. PG and LG peak amplitude arrival times were read automatically 
from segmented detector output files, and epicenter coordinates were obtained by a 
grid-search procedure (Ruud etal., 1993). Star is the mine site, diamonds, the envelope 
locations, and dots, the analyst-derived bulletin solutions. Almost all locations are 
within 10 km of the mining site, (b) Location of Geiranger quay construction explo¬ 
sions. Symbols as for (a). Station MOL (distance 55 km) is critical for accurate locations 
here, but only occasionally could PG arrival times be obtained due to LG interference. 
Also shown (white triangles) are new epicenter locations obtained by replacing PG 
times with corresponding Pg and Pn arrival times as read by analyst, (c) Rana mining 
explosions (two sites) used for preliminary envelope location experiments. Symbols as 
for (a)—same travel-time curves as for Titania. Since the network stations are deployed 
in the N-S direction (the ‘‘width'* of Norway at this point is only 5 km), the E~W 
location is not well constrained. Bulletin solutions often incorporated readings from 
one or more arrays to the east, giving better E-W resolution, (d) Svartisen dam exca¬ 
vation explosion used in envelope location experiments. Symbols as for (a). Note the 
excellent latitude estimates obtained by the envelope locations (diamonds). 
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Table 1 

Summary of location results for the four sites. Latitudes and longitudes are average results from a number of events. Standard 
deviations are computed relative to the average locations of the same events. Median location errors are computed relative to the known 
locations of the sites. Largest deviations are found for the N. Norway explosions Rana and Svartisen; the bulletin epicenter solutions 

benefit from access to phase readings from the arrays to the east. 


Site 

Method 

Latitude 
(deg. N) 

Std. Dev. 
laL (deg) 

Longitude 
(deg. E) 

Std. Dev. 
long.(deg) 

Median Loc. 

Err. (km) 

No. of 

Events 

Titania 

58.33 


6.43 




Bergen bull. 

58.27 

0.10 

6.41 

0.18 

11 

29 

Enve. P S 

58.43 

0.08 

6.31 

0.20 

17 

29 

Man. P + Enve. S 

58.26 

0.06 

6.40 

0.15 

9 

29 

Geiranger 

62.10 


7.21 



18 

Bergen, bull. 

62.12 

0.03 

7.29 

0.11 

6 

Enve. P + 5 

62.03 

0.05 

7.51 

0.20 

12 

18 

Man. P +• Enve. S 

62.09 

0.04 

7.29 

0.14 

8 

20 

Rana 

66.42 


14.69 




Bergen bull. 

66.45 

0.05 

14.70 

0.19 

11 

8 

Enve. P + 5 

66.33 

0.04 

15.12 

0.53 

15 

8 

Man. P -+ Enve. 5 

66.39 

0.06 

14.99 

0.20 

13 

16 

Svartisen 

66.77 


14.13 



12 

Bergen bull. 

66.81 

0.10 

14.50 

0.47 

. 8 

Enve. P + 5 

66.75 

0.02 

14.56 

0.92 

8 

12 

Man. P Enve. S 

66.72 

0.34 

14.58 

0.72 

15 

19 



Figure 5. (a) Record displays for a local earthquake in Bothnian Bay, Finland, and 

recorded by the center seismometer for several array stations (Fig. iX A 3 to 6 Hz 
bandpass filter was used, (b) The envelope transforms of the waveforms in (a) after 
smoothing, otherwise-as for Figure 2b_Note the.dominance. o£-LG_ancLthaLBG_is- 
reasonably prominent out to about 500 km. An alternative if PG is unclear is to read 
the first-arriving P wave (F/i). Our envelope epicenter solution coincided with that in 
the bulletin and shown in Figure 1. Phase readings from one or more of these arrays 
are used frequently in the Bergen bulletin epicenter solution; corresponding observa¬ 
tions for the other events in this study were not available for our envelope locations. 
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Iris Newsletter 


Comments on “The Political Sensitivity of Earthquake 
Locations by van der Vink and Wallace” 

M. Henger and K. Koch, Federal Institute for Geosciences and Natural Resources, Germany 
B. O. Ruud and E. S. Husebye, University of Bergen, Norway 


In the Fail 1996 IRIS Newsletter, van 
der Vink and Wallace (1996) use as one 
of their examples a small seismic event 
(magnitude 2.4) on 13 January 1996 
near the former USSR test site on 
Novaya Zemlya. Other seismologists 
who were requested to evaluate this 
event obtained a different epicenter 
location from that of van der Vink and 
Wallace although they used the same 
data, namely that from the ARCESS 
(ARAO) and Spitsbergen (SPAO) arrays. 
Data from the three-component Norilsk 
(NRI) station in Siberia indicates that 
the location determined by van der Vink 
and Wallace is incorrect and supports 
the other determinations (Figure 1). 


Data Analysis 

We extracted relevant time segments 
of recordings from the Norilsk three- 
component station. At ARAO and SPAO 
the noise conditions are dominated by 
strong, low-frequency microseisms up 
to 2 Hz, so the optimum signal-to-noise 
ratios (SNRs) are in the 2 -10 Hz band. 
At NRI the local noise conditions are 
very different with much less noise at ’ 
the low frequencies. The S-wave is 
visible even at very low frequencies and 
there is a clear wavelet near the expected 
P arrival time with a dominant frequency 
of 1.5 Hz (Figure 2a & b). 

The methods for analyzing array- 
recorded signals are well established 
and in cases with good SNRs even simple 
bandpass filtering and subsequent 
stacking suffice to accurately pick Pn 
and Sn arrival times. For the 13 January 
event, this was not the case, due to 
interference in the SPAO records with 
signals from an earthquake on the mid- 
oceanic ridge west of Spitsbergen 
(Ringdal, 1996). Van der Vink and 



Figure 1. Various epicenter determinations 
for the 13 January 1996 event at/near Novaya 
Zemiya: (1) two-array location (SPAO & 
ARAO) using lASPEI tables, (2) as for (1) but 
using Bergen University tables, (3) two arrays 
and the NRI station using lASPEI tables, and 
(4) same as (3) but using Bergen University 
tables. Note, RingdaPs (1996) two-array 
location at 75.38®N, 56.7®E is close to our 
epicenter (4). That of van der Vink and Wallace 
(1996) would coincide with (1) if we add 5 sec 
to their Sn-reading. The Novaya Zemlya test 
sites are marked by stars. The map insert 
shows the seismic arrays and 3C stations in the 
vicinity of Novalya Zemlya. Records from 
APAO (Apatity) and the new AMD (Amderma) 
array were not available to us. To our 
understanding the APAO array did not detect 
this event. 


Wallace, being unaware of this, picked 
a too early Sn arrival at SPAO, which 
resulted in a shorter epicenter distance 
and hence ‘moved’ the epicenter well 
into the Barents Sea. 

Apart from the misinterpretation of 
the Sn arrival due to the interfering 


event, another aspect contributed to the 
move of the epicenter westward. The 
differential epicenter angle between 
SPAO and ARAO is only about 45 deg, 
so the epicenter location resolution is 
relatively poor in the east/west direction 
as it depends strongly on the traveltime 
model used. However, the accuracy 
improves considerably if relevant Pn 
and Sn arrival times can be extracted 
from the NRI records. An additional 
advantage here is that with three station 
observations the choice of traveltime 
table used for location is less critical 
(Kennett, 1996) if the azimuth coverage 
is reasonable. 

Envelope processing (Husebye et al., 
1998) proved successful, as 
demonstrated in Figure 2c. The top trace 
is the envelope for the 1.0 - 2.5 Hz 
passband and a prominent, presumably 
P-wave arrives at 113 sec. It was also 
seen in the original records. However, 
extensive 3C analysis gave an azimuth 
of about 260 degrees which is far off the 
310 degrees expected for Novaya 
Zemlya and polarization characteristics 
were non-P. The corresponding 
waveform has some semblance to a local 
Rg-wave - no sharp onset and 
monochromatic pulses. In Figure 2c we 
also display the 3C envelope for the 2 - 
5 Hz passband; a presumed weak P- 
arrival at 122 sec (SNR -- 1.3) and a 
clear, presumed Sn-arrival at 243 sec 
(SNR ^ 1.5). A SNR of 1.5 is taken to 
indicate a significant signal in the 3C 
envelope since this process is similar to 
incoherent beamforming with a rather 
long STA window (a threshold of 1.6 is 
used for NORSAR incoherent beams). 
As an additional check, we used the 
Pisarenko et al. (1987) phase-picking 
algorithm, whose picks coincide with 
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Notch filtered and bondposs fiflered 1.0 - 2.5 Hz 



Notch filtered and bondposs filtered 2.0 - 5.0 Hz 



Figure 2. The 3-component NRI (Norilsk, 
Siberia) records for the 13 Jan 1996 event. The 
upper and middle parts) show filtered records 
in the 1.0 - 2.5 Hz and 2,0 - 5,0 Hz passbands 
respectively. Notch filters were used in addition 
to bandpass for removing spectral spikes. In c) 
the two upper traces are 3 C envelopes (Husebye 
et al. 1998) Pn-and Sn-arriva! times as picked 
for this event are marked. 

the envelope records (Figure 2c). Using 
the derived NRI Pn and Sn arrival times, 
we obtained an epicenter location of 
75.38° N, 56.55° E, using the Bergen 
University traveltime curves (Figure 1). 
The residuals using the Bergen 
University traveltime curves are not 
unreasonably large (less than 1.6 sec), 
while for the lASPEI model the residuals 
varied +/- 5.6 sec, although the 
corresponding epicenter difference was 
less than 15 km. Note, the weak Pn- 


arrival in the records are not critical for 
the epicenter determination; we have 
also located the event using a Pn arrival 
at 113 sec (see Figure 2c) and also 
without the P-arrival of NRI. In the first 
case the epicenter moved 10 km and in 
the second case only 2 km (in both cases 
the local traveltime tables were used). 

Discussion and Conclusion 

Reexamining the van der Vink and 
Wallace (1996) location, we found that 
their small misreading of the Sn-arrival 
time was not as critical as their use of 
the lASPEI tables, Non-Scandinavian 
scientists appear to consider the lASPEI 
tables to be adequate for signals recorded 
from Novalya Zemlya at stations in 
Fennoscandia, but the “local” scientists 
use somewhat different tables, which 
are more accurate for this region. 

What can be said about the source of 
this much debated event of 13 January 
1996? Its mb magnitude of only 2.4, 
which is equivalent to a yield not 
exceeding 100 tons TNT. Such small 
charges are typical of many chemical 
explosions and, even if it were nuclear 
in origin, it could hardly be termed a 
nuclear device in a CTBT context, as 
argued by van der Vink and Wallace, A 
criterion in favor of an earthquake source 
is the relatively small P/S-signal ratio 
(van der Vink and Wallace, 1996). We 
are not convinced that such a criterion 
can be used here without extensive 
modification, because the Barents Sea 
sedimentary basin blocks Lg- 
propagation, in other words, Lg-waves 
leak as S-waves into the upper mantle 
and reappear in the Sn coda (Mendi et 
al. 1997). Observationally, this is 
manifested in an extended Sn-wave 
train, as seen in Novaya Zemlya 
recordings. 

Using the novel envelope analysis 
technique we could extract useful 
traveltime information from the Norilsk 
3C station, despite the poor quality of 
the NRI records for an event of 
magnitude 2,4 at an epicenter distance 
of about 1200 km. However, we confess 


that such low magnitude events can 
neither be accurately located nor clearly 
identified with seismological means 
only. In particular, the LASPEI travel 
time tables are not well suited for 
accurate event locations in the Novaya 
Zemlya area for regional epicenter 
distances. Clearly, certain event location 
areas remain politically sensitive, as 
this event has demonstrated. 
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Response from van der Vink and Wallace 

Gregory van der Vink, IRIS, and Terry Wallace, University of Arizona 


We agree. Henger et al. provide a 
more accurate location for the epicenter 
of the 13 January event in Novaya 
Zemlya — one of the examples 
presented in our article. We based our 
location on the data from the ARCESS 
(ARAO) and Spitsbergen (SPAO) airays, 
and were unaware that signals from 
Norilsk (NRI) indicated a 
contemporaneous earthquake on the 
mid-ocean ridge west of Spitsbergen. 

The conclusions of our article, however, 
are in no way altered by the 
refinement of the epicenter 
location. As we have 
repeatedly seen, earthquakes 
that occur in sensitive areas 
during politically sensitive 
times are vulnerable to 
misinterpretation. 

Independent of which 
location is used (all of the 
epicenter determinations for 
the 13 January 1996 event 
are more than200kilometers 
from the Russian test site), 
the eaithquakeraisedconcem 
over Russian compliance 
with the ban on nuclear 
testing. Given the 
lowering of detection 
thresholds and the continuation of 
subcritical experiments, it is not 
surprising that such coincidences occur. 
And in fact, since the writing of our 
article, such a coincidence has occurred 
again, (see news articles inset and the 
bannergram of this Newsletter). 

To help prevent such false alarms in 
the future, data from all available stations 
needs to be examined. Henger et al.’s 
extraction of critical travel-time 
information for a magnitude 2.4 event 
from a single 3-component station at a 
distance of 1200 km is an enlightening 
example for those who may still debate 
the importance of auxiliary stations for 
the monitoring of the CTBT. For the 


more recent 16 August 1997 false alarm, 
critical data came from the IRIS GSN 
station KEV, which is not even part of 
the official CTBT monitoring system. 

We disagree with Henger at al in their 
assertion that a magnitude 2.4 event, 
which is equivalent to a yield not 
exceeding 100 tons of TNT, “could 
hardly be termed a nuclear device in a 
CTBT context”. From the United States ’ 
perspective, it most certainly would. On 
August 11, 1996 President Clinton 


rejected proposals for continued sub- 
kiloton tests, and announced his decision 
to pursue a “true zero-yield 
comprehensive test ban.” He based such 
a decision in recognition that the act of 
nuclear testing, not the threshold of such 
tests, is objectionable within the non¬ 
proliferation regime. 

Both the United States and Russia are 
conducting hydrodynamic experiments 
at their test sites. Whether these 
experiments must remain sub-critical, 
or whether they can produce modest 
yields is unclear from the negotiations. 
Because nuclear materials such as 
plutonium are used in the experiments, 
they are conducted in a manner similar 


to nuclear tests. Although there is 
currently no precise agreed upon 
technical definition of the maximum 
explosive energy from hydrodynamic 
experiments, a seismic signal 
corresponding to 100 tons of tamped 
explosive would be assumed as a 
violation of the ban on nuclear testing 
by the United States. 

Areas such as Novaya Zemlya, that are 
generally considered aseismic at 
teleseismic magnitude levels, wiH appear 
seismic at the lower 
magnitude levels detectable 
from regional coverage. 
Regional networks now 
provide a detection 
capability for the Novaya 
Zemlya area that is near mb 
2.5 (NORSAR, 1996). 
Accordingly, we can expect 
to find events of concern 
every few years, such as 16 
August 1997, 13 January 
1996,13 January 1995, and 
31 December 1992. Open 
networks in the western 
United States, where the 
seismicity is greater and the 
coverage more extensive, 
provide a detection 
threshold below magnitude 2.2 for the US 
nuclear test site (Hennet et al. 1996). In 
fact, six days before the US’ second 
1997 subcritical experiment, a seismic 
event was recorded on the Nevada Test 
Site. 
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SUMMARY 

The North Sea Lg blockage for wave paths across crustal graben structures is a well- 
established observational fact. Analysis of such observations implies that Lg blockage 
takes place in graben areas associated with sedimentary basin formation and crustal 
thinning. These intriguing observations have triggered many theoretical studies aimed 
at highlighting specific Lg loss mechanisms, albeit so far with only moderate success. 
Our approach to this problem is to simulate seismic wavefield propagation through the 
crustal waveguide using 2-D fimte-difference techmques. The graben structures are 
known in detail from oil exploration works in the North Sea, which has enabled us to 
use realistic crustal models in our Lg synthetics. In the most extreme model tested, the 
crystalline crust thickness beneath the graben amounted to only 5 km, while the over- 
lying sedimentary pile is nearly 10 km thick. At the base of the crust in the graben area 
the Moho is elevated nearly 10 km. This model has similarities to the oceanic crustal 
waveguide, where total Lg blockage is claimed for path lengths exceeding 100 km. The 
synthetic wavefields are displayed in terms of snapshots, semblance velocity analysis 
and time-space rms amplitudes. The dominant structural Lg loss mechanisms are the 
delay of the Lg waves in the thick sediments, Lg-io-Rg conversions (scattering) by 
lateral heterogeneities in the sediments, and 5-wave leakage out of the crustal wave¬ 
guide and into the upper mantle. A fraction of these upper-mantle S waves return to the 
crust and appear as Sn coda. Observationally, strong Sn phases of long duration 
are often associated with weak Lg phases and vice versa. Our synthetics produced 
Lg amplitude decay amounting at most to 6-10 dB, while observational data 
imply blockage amounting to 15-20 dB. The latter is equivalent to a Pn-Lg magnitude 
difference of nearly one magnitude urdt. The main outcome of this study is therefore that 
Lg-wave propagation is very robust and that a dominant blockage effect associated with 
intrinsic attenuation, that is Q values of the order of 100 at 2 Hz for a path length of 
minimum 100 km, is necessary to conform to observations. 

Key words: attenuation, finite-difference methods, Lg blockage, North Sea grabens, 
scattering, synthetic seismograms. 


INTRODUCTION 

An outstanding puzzle in observational seismology is that 
of Lg-wave blockage, which initially was reported for 
mixed ocean-continental paths (Press & Ewing 1952; Bath 
1954). More recently, Lg blockage has been reported for 
intercontinental wave paths, particularly across structural 
boundaries such as grabens, mountains and sedimentary 

_basimareas (Gregersen 1984;. Baumgardt 1990). In view of the 

importance of Lg observations for event magnitude and yield 
estimation, and also as a potential discriminant, many Lg 
studies have been published recently (Israelsson 1994; Zhang 
efa/.1994;Xie& Lay 1995). 

In view of the observational data at hand, many 
hypotheses or combinations thereof have been suggested 


to account for Lg blockage, namely crustal thickness 
variations, surface topography, geometry and thicknesses of 
sedimentary basins, mode conversions/leakage of Lg energy 
downwards into the mantle, Lg-to-Rg trapping in overlying 
sediments and naturally strong anelastic attenuation. In a 
recent study, Zhang et al. (1994) used a multivariate analysis 
approach in an attempt to quantify which of the above 
structural features is dominant for Lg blockage. Results were 
that waveguide geometry (abrupt crustal thickness variation) 
and anelastic attenuation were most important in this regard. 
Gregersen & Vaccari (1993) reached similar conclusions for 
Lg-blockage observations across the North Sea Viking 
Graben, and also pointed to the spreading of 5-wave energy 
over a longer duration due to interaction with the over- 
lying low-velocity sediments. Cao & Muirhead (1993) found 
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from finite-difference simulations that efficient Lg blockage 
occurred for models with a water depth of 2000 m. Their 
result does not apply to the North Sea grabens, since the water 
depth here (average ca. 100 m) is only a small fraction of the 
dominant wavelength. 

The anomalous £g-blockage observations have triggered 
many studies to explain this phenomenon in a theoretical 
manner. For example, Kennett (1986) used a ray diagram 
approach to explain qualitatively Lg blockage, while Bouchon 
(1981) used a discrete wavenumber representation for studying 
wave propagation in multi-layered media. Other scientists 
have studied Lg blockage effects across the Alpine Range 
and the Pyrenees (e.g. Chazalon et al 1993; Campillo et ai 
1993). Maupin (1989) used a coupled-mode approach for 
simulating North Sea Lg blockage, and Gregersen et al (1988) 
and Gregersen & Vaccari (1993) used a similar approach for 
analysing Lg blockage in the same area. In recent years, finite- 
difference and finite-element methods have become popular 
in computing realistic wavefield responses as a tool for 
explaining Tg^-blockage mechanisms, and we refer here to 
Regan & Harkrider (1989), Cao & Muirhead (1993) and Xie 
& Lay (1994). The feature in common for all the theoretical 
studies mentioned is the inability to explain fully the Lg 
blockage as observed across grabens in the North Sea and in 
other areas using realistic crustal models. The only way to 
obtain sufficiently strong blockage seems to be to introduce 
strong anelastic attenuation, since deformation of the crustal 
waveguides is unable to do this. As mentioned, for mixed 
continental-oceanic wave paths where the length of the latter 
is 100 km or more, Lg blockage is almost complete. Recently, 
Zhang & Lay (1995), using a normal-mode Lg representation, 
claimed that the overall thickness of the oceanic crustal wave¬ 
guide is too thin to allow Lg to develop as a significant phase 
in the frequency band 0.3-2.0 Hz, thus explaining oceanic Lg 
blockage. 

A general drawback with most of the theoretical approaches 
mentioned above is that the crustal models used are lacking 
many details on Moho undulation and sedimentary basins 
along the Lg propagation paths. Naturally, finite-difference 
modelling can easily incorporate such structural features 
and thereby provide better insight into Lg loss mechanisms. 
There have already been some efforts along these lines, such 
as the work of Cao & Muirhead (1993), Xie & Lay (1994) and 
Regan & Harkrider (1989) demonstrating in a convincing 
manner the complexity of Lg propagation even for slowly 
varying crustal models. In this paper, we simulate Lg blockage 
across the Central Graben through 2-D finite-difference 
(FD) synthesis, since Lg observations are readily available 
(Gregersen 1984; Kennett & Mykkeltveit, 1984; Kennett 
et al 1985) and because the graben structure is mapped in 
considerable detail (Barton & Wood 1984; Gabrielsen et al 
1990). 

OBSERVATIONAL DATA 

The Lg blockage here is well documented in the literature, 
in particular in the works of Gregersen (1984), Kennett 
& Mykkeltveit (1984) and Kennett et al (1985). The most 
pronounced blockage takes place for Lg paths with source 
locations to the west of the graben and landward recordings 
in SW Norway. To illustrate the North Sea Lg blockage, a 
few events have been chosen (source and receiver geometries 
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Figure 1. The two dominant North Sea graben segments are outlined; 
Central Graben to the South and the Viking Graben to the North. The 
seismograms displayed in Figs 2-4 are recordings from the earthquake 
locations marked by open circles. Station locations are marked with 
filled circles on the map. 


in Fig. 1), and the seismograms displayed in Figs 2-4. The 
graben blockage effects are obvious; outstanding features are a 
relatively strong Sn phase and a nearly extinct Lg phase. The 
expected Lg arrivals are literally hidden in the Sn coda or tail in 
most of the seismograms. 

MODEL BUILDING STRATEGY 

The amplitude patterns are not always consistent, as small 
source-receiver path changes could exhibit large differences 
in Lg propagation efficiencies (Mykkeltveit & Husebye, 
1981). Such variations are also obvious from North Sea 
event recordings as displayed in Figs 2-4. This observational 
variation in Lg blockage, part of which may relate to the source 
mechanism and focal depth, motivated a flexible modelling 
scheme as illustrated in Fig. 5. The idea here is that by sub¬ 
dividing the crust near the free surface and close to the Moho 
we can easily change the model types from the extreme of one 
layer over a half-space to others with thick unevenly distri¬ 
buted sediment strata, Moho ramps and associated crustal 
thinning. This approach allows great flexibility in model 
building; details are given in Table 1, which lists the set of 
graben models used in this study. Model 2 is the proper ‘graben 
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Figure 2. Seismic recordings from an earthquake in SE England 
(marked as 1 in Fig. 1) on 1994 February 15. The stations (see Fig. 1 for 
location) are operated in an automatic trigger mode, and the first- 
arriving P waves are occasionally lost. The expected arrival times of the 
Pn, Pg, Sn and Lg phases, calculated for the hj^pocentres and earth 
model used in the Norwegian Seismic Network Bulletin, are marked 
on each seismogram (for Lg a constant group velocity of 3.56 km s“^ 
is used). The traces have been bandpass filtered between 2 and 4 Hz. 
Outstanding features are strong Sn phases with long coda. The absence 
of the Lg phase is presumably caused by the North Sea graben block¬ 
age. The magnitude (M^) of this event is reported as 4.0 by British 
Geological Survey. The Pn and Lg magnitudes for the displayed 
seismograms are estimated to be 3.9 and 2.7 respectively. 


Figure 3. Seismic recordings from a central North Sea event on 1995 
June 28 (marked as 2 in Fig. 1). Filtering and phase markings as for 
Fig. 2. Outstanding seismogram features for the first five stations 
(from EGD to HYA) are clear Sn phases of moderate duration 
followed by more prominent Lg phases. The far-away stations MOL 
and NSS are anomalous in this regard; strong Sn phases of long 
durations followed by relatively weak Lg wavetrains. The graben 
blockage appears modest for this earthquake, while some extra¬ 
ordinary structures between HYA and MOL/NSS cause Lg blockage 
at the latter two stations. Ml is reported as 3.2 by NORSAR. For the 
displayed: seismograms, Pn and Lg magnitudes are estimated to be 3.0 
and 3.1 except for the distant stations MOL and NSS which have Lg 
magnitudes of 2.4. 


model’, while Model 1 is a reference model with normal crust 
used for comparison (see also Fig. 10). Model 4 is exceptional 
as the crystalline crust thickness is reduced to about 5 km 
between the Moho bump and overlying graben sedimentary 
infill, and might thus be considered typical of oceanic crust. 
Model 3 is a reference model for Model 4. In Model 5 we have 
introduced random perturbations of velocities and density 
through von Karman functions (order 0.3) with rms velocity 
fluctuations of 8 per cent and correlation distances (horizontal 
and vertical) of 10 and 2.5 km (details in Hestholm et al 1994) 
in the central part of the model (see Fig. 'll). Apart from" this, 
Model 5 is similar to Model 2. The ‘prototype’ model (Fig. 5) is 
based on a cross-section of the Viking Graben published by 
Gabrielsen et al 1990, 


FINITE-DIFFERENCE SYNTHESIS 

Our approach to 2-D FD modelling, including free-surface 
boundary conditions and source representation, is detailed in 
Hestholm et al (1994). Absorbing boundary conditions were 
included in the form of 10 km thick damping layers along the 
side wails and bottom of the model. This technique ejBfectively 
prevents artificial reflections from the model walls. Due to the 
absorbing boundaries, the useful distance range is 480 km 
(instead of the original 500 km shown in Fig. 5) and for this 
interval the vertical component of the wavefield is extracted at 
1 km equidistant points along the top of the model. The source 
is located at a depth of 10 km and produces signal energy in the 
0-5 Hz band with a peak frequency of 2.5 Hz. We have used a 
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Figure 4. Seismic recordings from a northern North Sea earth¬ 
quake on 1995 April 21 (marked as 3 in Fig. 1). Filtering and phase 
markings as for Fig. 2, Outstanding seismogram features are strong 
Sn phases of long duration followed by weak Lg phases. Apparently, 
the graben blockage causes Lg wave energy to leak as S waves into 
the upper mantle and subsequently reappear in the crust as part of the 
Sn coda. Pn and Lg magnitudes are measured to be 2.0 and 1.6, 
respectively. 


rotation source that generates only S waves. However, F waves 
created by P-to-S conversions are also seen in the synthetics. A 
problem in many Lg-blockage studies is that of proper ampli¬ 
tude reference levels. Since the source is identical in all our 
synthetics, we can directly compare amplitudes for the various 
models. Since structural differences may give slightly different 
arrival times of the phases, the reference level could be slightly 
distorted. In practice, this is not much of a problem since the 
abundance of synthetic traces allows us to smooth the rms 
amplitudes over time and distance. In an analysis of real 
observations, the problem is more severe and an additional 
disadvantage is that potential reference phases such as Pn, Sn 
or even Pg themselves are distorted due to differences in source 
and/or propagation effects.,The amplitudes of these phases 
depend on non-stationary parameters such as path, source 
type and focal depth and are thus not suited as reference 
levels. Our 2-D FD scheme for the synthetic wavefield calcu¬ 
lation does not incorporate intrinsic attenuation (Q) and is in 


Distance [km] 



Figure 5, Prototype of northern North Sea graben model (Viking 
Graben). At the top there is a 2 km thick sediment layer extending over 
the entire length of the model. For the graben area (in the distance 
range 160-420 km) we have options for including the outlined sedi¬ 
mentary basins. Likewise, the blocks just above the Moho may be 
assigned either upper-mantle or crystalUne-crust physical properties. 
Further detaUs of the models actually used are given in Table 1. The 
source (open circle) is located 10 km from the left edge of the model 
and at a depth of 10 km. Note the different scaling of horizontal and 
vertical axes in this plot. In FD computations the model was extended 
downwards to a depth of 80 km. 


this respect similar to many other blockage studies reported in 
the literature. However, intrinsic attenuation can be estimated 
in a manner similar to that used in Lg magnitude studies and 
will be discussed in a later section. 

MODELLING RESULTS 

Visual inspections of seismic records are instructive and 
provide us with a rough insight into wavefield propagation 
complexities: However, with finite-difference modelling the 
wavefield can also be displayed as snapshots and the syntheUc 
seismograms can be analysed in terms of apparent velocities 
(semblance analysis) and rms amplitude variations. 


Table 1. Summary of models used in this paper. The prototype 
model, shown in Fig. 5, was manipulated by assigning specific sedi¬ 
ment, crust or upper-mantle physical properties to the various blocks. 
The F velocities for sediments, crust and upper mantle are 3.6, 6.7 and 
8.2 km s“^ respectively. S velocities are derived from P velocities 
using a Poisson’s ratio of 0.25, while densities are derived from P 
velocities using Birch’s law. ‘Crustal Thinning’ implies that the middle 
part of the crust is made 9 km thinner than in the ‘prototype’ model 
while otherwise retaining the model geometry. For models marked 
with ‘Sediment Layers’ all the blocks in the upper part of the model 
were given sediment properties. The 2 km thick sediment layer at the 
top was retained in all models. ‘Moho Bump’ means that the three 
blocks near the Moho were assigned upper-mantle physical properties. 
‘Velocity Perturbation’ of 8 per cent rms in the graben area was used 
only for Model 5; see Fig. 11 and text. 


Model 

Crustal 

Sediment 

Moho 

Velocity 


Thinning 
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Bump 

Perturbation 
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Seismograms 

In contrast to real seismograms (Figs 2^), even for the 
most extreme graben structures used in this study, the wave 
energy in the expected phase-velocity window for Lg, i.e. 
3.2-3.9 km s“^ is stronger than that for Sn. In the Hestholm 
et al (1994) study, we found that both Pn and Sn were always 
relatively weak phases, as expected from the strong geo¬ 
metrical spreading effects of head waves, being roughly pro¬ 
portional to where r is the epicentral distance. Somewhat 
stronger Pn and Sn waves were obtained by introducing 
velocity gradients and heterogeneities beneath the Moho, so 
that Pn and Sn become a mixture of head waves and scattered 
waves. For ray paths entirely within the Baltic Shield, the Sn 
phase is generally weak and often not reported by the analyst 
(Vogfjord & Langston 1990). 

A comparison of real (Figs 2-4) and synthetic records 
(Fig. 6) shows that observed Sn is relatively strong and of long 
duration, clearly implying that a considerable amount of shear 
waves leak out of the crustal waveguide (Gregersen & Vaccari 
1993). In contrast, the synthetics exhibit a relatively weak Sn 
and relatively strong Lg of long duration. This does not mean 
that there is no Lg blockage in our synthetics. First, synthetic 
sensor positions may be too close to the graben for leaking Sn 
phases to be observed. In addition, there is a considerable 
contribution of Rg waves, generated by scattering in the sedi¬ 
mentary basin, to the Lg wavetrain in our synthetics. Note that 
Rg waves are generally not observed on stations in Norway for 
paths longer than 100 km. The reason for this is that near¬ 
surface topography and heterogeneities, in combination with 
strong intrinsic attenuation in the uppermost 2 km of the crust, 
eflEiciently wipe out Rg waves for longer paths (Ruud et al 
1993). 

Snapshots 

Such plots, showing the total spatial wavefield for a specific 
time, are convenient for visualizing scattering phenomena, 
and, in our context, also the leaking of energy out of the crustal 
waveguide. The Model 2 synthetic wavefields are instructive for 
illustrating these phenomena, as shown in Fig. 7. The strongest 
amplitudes in these snapshots are associated with the crustal S 
waves, i.e. Sg and Lg. Further comments for each time frame 
are given below. 


at the surface. This result is due to the combined effects of 
the incoming SV, the reflected SV and the inhomogeneous 
P waves (see Aki & Richards 1980, p. 190, Problem 5.6). 
Although our model has a thin sedimentary layer near the 
surface, this result still seems to be approximately valid. The 
largest amplitudes usually recorded in the S wavetrain, that is 
the Lg phase, are created by S waves that have been repeatedly 
reflected within the crust (‘guided waves’)- These waves have 
higher apparent velocities but lower group velocities than the 
Sg phase. The waves with the highest apparent velocities, still 
undergoing total reflection at the crust-mantle boundary, will 
have the same apparent velocity as the Sn phase, i.e. about 
4.7 km s“^ in our model. These waves will terminate the Lg 
wavetrain since S waves with higher apparent velocities will 
loose part of their energy when reflected at the Moho interface. 
The ‘group’ velocity (distance divided by traveltime) of the 
latest Lg waves can be found as follows. The sine of the 
incidence angle is 3.9/4.7=0.82, The horizontal component of 
velocity of such a ray will thus be 3.9x0.82 = 3,2 km s“Mnthe 
analysis of Lg amplitudes we have taken this group velocity 
to define the end of the Lg window. The upper-mantle S wave 
(the Sn phase) has propagated out to about 270 km at this first 
snapshot and travels faster than the crustal S waves, reaching 
the edge of the displayed section at the 100 s time frame, as can 
be seen in the subsequent time frames. 

Time frame 70 s 

The squeezing of the crustal S waves associated with thinning 
of the crystalline crust is obvious in this time frame as 
well. Some energy is transferred into the overlying sediments 
(220-260 km) and in the later part of the Lg wavetrain much of 
the energy leaks down into the mantle. The dipping Moho 
interface results in a smaller incidence angle of the rays here so 
that the S waves are only partially reflected (instead of total 
reflection which is the rule for guided waves such as Lg). 

Time frame 80 s 

The direct Sg waves have now passed the thinnest part of the 
crystalline crust and relatively little Lg wave energy remains 
in the crustal waveguide, while forward-propagating S waves 
under the uplifted Moho are relatively strong. In the sedimentary 
basin much of the Sg and Lgv/^yQ energy seems to be trapped. 


Timeframe 60 s 

In the first snapshot, the P phases are seen in the distance 
interval 390 to 440 km. (The Pn appears stronger than Pg here 
because the latter have mainly horizontal particle motion—not 
shown in the snapshots) A strong Lg wavetrain has developed 
and has just reached the graben structure. The direct S wave 
{Sg) propagates at a velocity of 3.9 km s“^ and is seen at about 
235 km offset. The Lg wavetrain, observed as up- and down¬ 
going wavefronts, is spread over an interval from about 190 km 
and up to the Sg phase. While the Sg phase is the strongest 
one seen in the snapshots, it is rather weak in the synthetic 
seismograms .extracted at the free surface (Fig. 6). Taking a 
closer look at the snapshots, it appears that the Sg phase is 
hardly ‘reaching up’ to the surface. In fact, this is consistent 
with the theory that 2 in SV wave propagating parallel to the 
free surface of a half-space should result in zero amplitude 


Time frame 90 s 

Direct S and Lg waves are still confined to the thinnest part of 
the crustal waveguide (see Fig. 5). Between 350 and 380 km the 
crystalline crust thickens again and some of the sub-Moho 
propagating shear waves leak back into the crust. Rg waves are 
particularly energetic in the sedimentary basin, having phase 
velocities between 2.0 and 2.5 km s“^ 

Time frame 100 s 

The Sg wavefront has partly been restored and again fills most 
of the crustal waveguide. Some of the S waves, which h^aye 
leaked down into the upper mantle due to crustal thinning, also 
propagate approximately horizontally and partially leak up 
again into the crust as Sn coda waves (see the rather strong 
upper-mantle S waves around 400 km). 
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Model 1 Model 2 


20 40 60 80 100120 140 160 180 20 40 60 80 100 120140 160 180 



Figure 6. Synthetic seismograms for models 1 and 2 for the distance range 140-480 km using a station spacing of 20 km. The two dotted lines 
indicate the Lg group velocity window of 3.9 and 3.2 km s“^ The numbers on the left are the distances from the source in kilometres. Up to about 
150 km the two models are identical and the seismograms are very similar except for the backscattered waves seen in the later part of the Model 2 
records. From 160 and up to 340 km the amplitude level in the Lg window is larger in the graben model (2) than in the reference model (1). There are 
several reasons for this: (i) the increased sediment thickness tends to amplify the waves, (ii) some body waves are converted to Rg waves at lateral 
heterogeneities near the surface, and (iii) some waves that are squeezed out of the crustal waveguide continue to propagate in the sediments as guided 
waves. Beyond 340 km, the crust gradually regains its normal structure and the partial blocking of the Lg phase becomes evident. The wave energy 
which was forced up into the sediment layers is delayed relative to the Lg phase and becomes smeared over a longer time interval. In both models the 
Sn phase is weak compared to Lg. 


Time frame 110 s 

For this timestep we show the wavefields for both Model 2 and 
Model 1 (lowest section). Note the different outlines of the 
Moho interface and sedimentary basin for the two sections. 


The Sg and Lg waves of Model 2 have no w left the graben 
structure and some blocking effect is clearly seen when com¬ 
paring the two snapshots. For Model 1 the S wavefield is 
quite similar to that seen in the first snapshot, except for a 
horizontal stretching effect. The later part of the *S'-wavetrain 
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{Lg) has suffered more blockage than the Sg phase, which has 
propagated nearly horizontally in the middle part of the 
crust. By measuring amplitudes at 20 km depth, we find about 
6 dB blocking for Sg and 9 dB for Lg. At the free surface the 
blocking is not so obvious; the Sg amplitude is very low for 
both models and the Lg waves are partially hidden by the 
strong surface waves confined to the sedimentary layer. 

The essence of this series of snapshots is that when the S 
waves are ‘forced' into the graben thinning areas, part of their 
energy is ‘squeezed' out of the crystalline crust through S-to-Rg 
conversions in the sediments, and by shear waves penetrating 
the Moho interface. With increasing time and after the graben 
thinning area is overtaken, some S waves may leak back into 
the waveguide. In other words, a sort of Lg wavefield healing 
process appears to commence. 

Semblance analysis 

The synthetic records were also subjected to semblance 
analysis in specific distance intervals. This is essentially a 
time-velocity analysis convenient for the identification of 
various phases in the seismic records based on waveform 
similarities. It is computed as the normalized ratio of coherent 
energy to total energy within short moving time windows for 
varying stacking velocities (Taner & Koehler 1969). Semblance 
results for Model 1 (smooth) and Model 4 (complex) are dis¬ 
played in Figs 8 and 9 for two distance ranges: 240-250 km 
and 300-310 km. In the first case, the smooth crustal wave¬ 
guide model does not distort the synthetic wavefield much, as 
seen from the high semblance values signifying very good 
spatial correlation. Pn and Pg phase velocities in the velocity 
window 6.7-8.1 km s”^ are clearly observable in the time 
interval of 40-60 s with a subsequent Sn arrival at 62-67 s. 
These phases are followed by prominent Sg and Lg wavelets 
throughout the 67-90 s interval. The complex Model 4 
wipes out signal correlations both for P and S, although 
we see some P and S scattered phases in the interval 42-52 s 
and S-lo-Rg and P-to-i?^ scattered phases (velocity range 
2.0-2.5 km s“^) in the interval 52-90 s. Semblance results in 
the distance range 300-310 km (Fig. 9) are quite similar to 
those at 240-250 km (Fig. 8). The Model 1 results are essen¬ 
tially unchanged except for a 10 sec time shift while those of 
the complex Model 4 show more 5-type and Pg-type scattering 
contributions. 

rms amplitudes 

An instructive way to demonstrate how the seismic energy 
propagates across the model is given in Fig. 10 (Model 1— 
reference) and Fig. 11 (Model 5—^perturbed graben model). 
Here the rms amplitude (in dB) for a 2 s moving time window 
is displayed as a function of distance and traveltime. Pn 
energy is reasonably well preserved in both models, while the 
distorted crustal waveguide beneath the graben (Model 5) 
effectively prevents Pg-type propagation (Fig. 11). Sn phases 
are visible in both models, and direct Sg waves are clear out 
to about 325 km. Lg energy is seen in the group-velocity 
interval 3.2-3.9 km s”^, while in the interval 2.0-3.0 km s“^ 
•S-to-Pg scattering contributions dominate. After passing over 
the graben (beyond 400 km), it is obvious that some of the 
Lg energy has been blocked for Model 5 as compared to 
Model 1. There is also a considerable amount of wave energy in 
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the upper-left part of Fig. 11 (but not in Fig. 10), which is 
identified as backscattered Pg waves. 

In all discussion of Lg propagation, the term blockage 
is seldom given a quantitative definition. As mentioned, in 
the observational studies Pn and/or Sn phase amplitudes are 
used as a reference for measuring the extent of blockage. As 
demonstrated above, both Sn and in particular Pg phase 
amplitudes are affected by wavefield propagation through a 
complex waveguide. Our measure of relative Lg amplitude 
decay comes from comparing the root-mean-square (rms) 
amplitude distribution derived from the models with pro¬ 
nounced graben structure (Table 1) with those obtained for the 
reference model. The results are displayed in Figs 12 and 13 and 
the following comments apply: Pn amplitude decay is only a 
few decibels, while Pg amplitudes are down more than 12 dB 
or a factor of 4. Sn is down by a few to 6 dB, while the Lg 
waves, velocity interval 3.2-3.9 km s”^ are down 6-9 dB, i.e. 
a reduction by a factor of maximum 3, which is less than that 
for the Pg phase. In Figs 12 and 13, the slow Pg scattering 
wavelets appear as strongly amplified. The reason for this is 
that the reference models are almost void of Pg scattering 
contributions. 

It was a bit disappointing to us that, even with reahstic 
but nevertheless complex North Sea graben models, we are 
not able to synthesize Lg blockage as implied by real obser¬ 
vations. The Lg amplitude reduction obtained is similar to that 
reported by Zhang & Lay (1995) in their normal model simu¬ 
lation of oceanic Lg blockage which they deemed successful. 
For our part, we do not entirely agree with this, since an 
amplitude reduction factor of 3 to 4 is hardly equivalent to 
blockage. In this respect we tend to agree with Maupin (1989), 
who also failed to simulate Lg blockage satisfactorily for 
the North Sea structure and thus concluded that Lg wave 
propagation was indeed robust. A potential defect of our 
waveguide models is that wave velocities are homogeneous 
within individual sedimentary layers, with constant velocities 
of respectively 6.7 and 8.1 km s”^ in the crystaUine crust and 
below the Moho. As demonstrated by Hesthohn et al (1994), 
coda waves observed in local explosion and earthquake records 
require small-scale rms velocity perturbations of the order of 
2-4 per cent. Since graben structures exhibit more structural 
variations than crustal layers, we introduced rms velocity 
fluctuations of 8 per cent (see Figs 11 and 13) to explore to what 
extent this would affect Lg propagation. The outcome is that 
Lg blockage amounted to roughly the same as for a very thin 
graben crustal waveguide, as implied by a visual comparison 
between Figs 12 and 13. A summary for this result section is 
that Lg blockage for realistic North Sea graben strucUires 
seldom exceeds 10 dB, even for the most extreme models. In 
other words, Lg propagation is very robust even under adverse 
structural conditions typical of graben areas. 

In the above discussion we have identified various causes of 
jLg-blockage, such as (i) energy escaping through the Moho, 
(ii) energy delayed or converted in the sediments, or (iii) the 
thinning of the crystalline crust per se. It would be interesting 
to assess the relative importance of these mechanisms. 
Obviously, they cannot be considered to be independent of 
each other - as, for instance, the thinning 'of the" crystalline^ 
crust is a consequence of the thick sediments and the Moho 
shallowing. In physically realistic models, thick sediment 
basins are usually accompanied by thin crust (for isostatic 
compensation), but in our modelling experiments we can 
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Figure 8. Semblance analysis results for the distance range 240-250 km that is the graben entrance (see Fig. 10). The analysis was performed with 
bandpass-filtered (0.5-2 Hz) traces and a 1.5 s moving time window. The amplitude shown above each semblance plot is relative to the maximum 
rms amplitude of each section. For Model 1 (reference model—normal crust) the wavefield distortions are modest so the signal correlations across 
the 11 traces in this 10 km line array are good. Until the Lg waves start at about 65 s the synthetics are dominated by waves with apparent velocities 
above 6 km s“^, i.e. F-type velocities. As expected, the apparent velocity in the Lg wavetrain is seen to increase slightly with time from about 4.0 to 
4.2 km s*^^ For Model 4 (graben model with thin crust) the wavefield is severely distorted so semblance is weak. For this model the scattering contri¬ 
butions in terms of S (velocities about 4 km s “ ^) and Rg waves (phase velocities slightly above 2.0 km s~ ^) are seen over almost the entire time interval. 


ignore this principle. Thus, besides the models in Table 1, we 
have tested a model with thick sediments but without Moho 
shallowing and another model with a shallow Moho but 
without thick sediments. These experiments indicate that thick 


sediments are a more effective blocking mechanism than 
crustal thinning. More interestingly, however, the sum of the 
individual blocking mechanisms (as measured in dB) is less 
than when they are working together. This result is reasonable: 
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Figure 9. Display of the semblance analysis results for the distance range 300—310 km, which is the central part of the graben structure for models 1 
and 4. Results are quite similar to those in Fig. 8. 


when the crystalline part of the crust gets thinner the Lg waves 
must undergo more reflections to traverse the same horizontal 
distance. Every time the waves are reflected from the free 
surface^ they must pass through the sediments-(except.for the 
horizontally propagating Sg wave which is reflected from the 
bottom of the sediments). The delay these waves experience 
depends mainly on the sediment thickness and number of 
reflections. The conversion of Lg waves to I^g (or other waves) 


in our model is probably controlled by the ‘roughness’ of the 
crust/sediment interface. We have not varied this parameter, so 
we cannot judge its importance. However, our model appears 
rather rough in this respect so it is unlikely that this effect is 
underestimated in our modelling. Lg wave energy ‘penetrating’ 
the Moho is most likely to occur for S waves with high 
apparent velocities, that is, for the later part of the Lg wave- 
train. For a Moho incidence angle of 10° less than the critical 
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angle, about 90 per cent of the energy will escape into the 
mantle (computed from plane-wave energy transmission/ 
reflection coefficient). However, this mechanism is effective 
only while the Lg waves are propagating ‘updip’, whereas the 
effect of the thick sediments is active over the complete range of 
the graben structure. 

To summarize, we find that the combination of thick 
sediments and thin crust is the most effective structural Lg 
blocking mechanism. This is confirmed by another result of 
our FD modelling experiments: if we compare the rms'ampli¬ 
tudes of Model 2 with those of Model 1, we find that the Lg 
blockage is about 3 dB less than when comparing Model 4 
with Model 3 (Fig. 12). The only difference between these 
two cases is the thickness of the crystalline part of the crust; 
the thicknesses of the sediments and the Moho structure 
are exactly the same. We also note that Model 2 is the one 
which best fits the North Sea grabens with respect to layer 
thicknesses. 

DISCUSSION 

From the results presented in the previous section, it is clear 
that even graben models with more adverse crustal thinning 
and velocity variations than what is found in the North 
Sea area are unable to block Lg waves to the extent observed. 
This also applies to cases where Lg traverse crystalline crust 
as thin as oceanic crust. With respect to model complexities, 
our modelling approach is very flexible in the sense that 
crucial crustal waveguide parameters, i.e. thickness, extent 
of sediment structures and Moho roughness, can easily be 
modified. In other studies, the Moho is often described by 
a simple ramp function and the sediment infill above by 
a trough. As mentioned, Lg blockage was not very efficient 
for our models—a factor of 2-to-4 in Lg amplitude reduction— 
so we introduced small-scale structural complexities in terms 
of 8 per cent rms velocity perturbation in the crust beneath 
the graben. The additional blockage effect of such extreme 
structures amounts to a maximum of 3 dB. In other words, 
Lg blockage cannot be fully attributed to structural features 
alone so the effect of intrinsic attenuation (Q) must be of 
importance. As mentioned, it is not a trivial task to incorporate 
such effects in a FD modelling scheme, although it can be done 
(Robertsson et al. 1995). However, since intrinsic attenuation 
is essentially a path-length effect, we can estimate fairly 
accurately the Q necessary to give a total Lg blockage of 
20 dB relative to a ‘normal* crust. In this regard we attribute 
10 dB Lg amplitude loss to complex graben structures and 
need to account for an additional 10 dB amplitude reduction 
over a graben width of 100 km. In this regard, a Q value of 
50 at 1 Hz (or 100 at 2 Hz) is needed for a 10 dB amplitude 
reduction. In continental crust, Q values are much higher 
but such values are not extreme in graben areas that have 
been subjected to much deformation in terms of stretching, 
thinning and block faulting. Literally a ‘joker* in many geo- 
dynamical problems is the extent of water within the crust/ 
upper mantle, since a water content of just a few per cent may 
significantly change physical rock cha^racteristics (Anderson 
-1990). In the case of attenuation properties, rock inclusion of 
even moderate water content would tend to lower the Q value. 

Of course, it is not unreasonable to expect that crystalline 
rocks in oceanic and marine shelf areas have a higher water 
content than in continental areas. In other words, a 10 dB 


amplitude reduction of Lg waves propagating through oceanic 
and graben areas of 100 km extent appears quite realistic. 

At this stage, it may be appropriate to return to real North 
Sea seismic event recordings; a few examples are given in 
Figs 2-4. A remarkable feature for events 1 and 3 is the strong 
Sn phase, which, in addition, often has a long duration. Similar 
observations were reported by Baumgardt (1990) for events in 
Novaya Zemlya with paths crossing the Barents Sea. In con¬ 
trast, all Canobe profiling shots in the North Sea (Fig. 2 
m Kennett & Mykkeltveit 1984) exhibit weak Sn phases, 
irrespective of shot location west or east of the graben 
structures. For the three shot points west of the graben, the Lg 
phase is hardly visible compared with shot points east of the 
graben. The recordings here were at the NORSAR array to the 
east of the graben. Returning to our observations. Event 2 
(Fig. 3) exhibits some remarkable features for an epicentre on 
the western side of the graben, namely that Lg is strong for the 
nearest stations on the west coast of Norway, while for stations 
further north (MOL and NSS) Sn is relatively strong and Lg 
correspondingly weak. This last point is the more surprising 
since the source-receiver paths to stations SUE and HYA 
initially coincide with those to MOL and NSS, so partial 
Lg blocking apparently takes place along the on-land path 
segments where sediments are lacking. For the MOL and NSS 
stations we notice the observational feature typical of events 1 
and 3, namely a relatively strong Sn phase of considerable 
duration. These Sn features are typical for cases where Sn 
is relatively strong and Lg correspondingly weak, so it is 
tempting to explain this in terms of Lg leakage from the crustal 
waveguide into the mantle where subsequently some of it is 
returned to the crust and appears at the free surface as delayed 
Sn phases. The mode of returning 5'rwave energy leakage 
from the crustal waveguide is demonstrated in our snapshots of 
the synthetic wavefield (Fig. 7). The example with Sn~Lg 
recordings for stations MOL and NSS for Event 2 clearly rules 
out the strong Sn phase being a specific source effect. In 
general, for most local event recordings in Fennoscandia the 
Sn phase is at best weak and sometimes not observable at all. 

In this respect, the strong Sn phases often, seen in North Sea 
event recordings could be a result both of strong attenuation 
of the Lg waves in the crust and of excessive ^'-wave energy in 
the upper mantle due to Lg leakage. We may have simulated 
this phenomenon more realistically by introducing a strong 
velocity gradient just below the Moho so that Sn becomes a 
combination of direct (head) and diving waves and thus 
appears stronger (Hestholm et al 1994). 

Finally, we also examined some of the observational 
evidence reported in the literature for Lg blockage in the North 
Sea. When Lg amplitude decay is estimated from the Sn/Lg 
ratio there may be a bias since Lg energy leaking into the 
mantle seemingly returns later in the Sn coda. Another aspect 
here is that in the North Sea Lg blockage and graben structural 
complexities are almost synonymous concepts, which in turn 
reflect observational results presented by Gregersen (1984) and 
Kennett et al (1985) among others. Kennett et al (1985) used a 
back projection scheme for identifying blockage areas in the 
North Sea and found these to be coincident with the dominant 
Central and Viking grabens (Fig. 1), A visual inspection of the 
observational data at hand gave that the wave path cross¬ 
sampling was poor; that is Lg-blockage contributions may also 
reflect off-graben structural peculiarities. In essence, graben 
blockage of Lg waves may be less than generally assumed and 
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Figure 7. A series of snapshots illustrating wave propagation through a typical graben area (Model 2 in Table 1). The horizontal axis gives the 
distance from the source (in km) and the vertical axis gives the depth from the surface (in km). The Moho and the sediment/crust interfaces’are shown 
as solid lines. The time (in s) is given in the lower left comer of each snapshot. The two last snapshots are for the same time but for different models in 
the lowermost snapshot, the wavefield of Model 1 is shown for comparison. See text for discussion. 
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Figure 10. RMS amplitude display for Model 1 (the reference model). The lower part of the figure shows a section of the model. RMS ampli¬ 
tudes were computed for each sensor at the surface (1 km spacing) with a 2 s long moving time window. On the right vertical axis of the rms plot 
the group velocities are given in km s”^. T 3 ^ical Lg group velocities are in the interval 3.2 to 3.9 km s“^ and coincide with the strongest amplitude 
in this plot. 
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Figure 11. RMS amplitude display for Model 5 with 8 per cent velocity perturbation of the crustal velocities in the graben region as shown in the 
lower part of the figure. Pg and Lg amplitudes are reduced compared to Model 1 (Fig. 10), while the outstanding features are the large amplitudes 
stemming from forward and backward S-to-Pg scattering from the sedimentary basins and in particular from the perturbed region. The group 
velocity of these Rg waves (propagating mainly in the sedimentary layer) is less than 2.0 km s“^ 
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Figure 13. Relative rms amplitudes of the most complex crustal model (5, Fig. 11) considered in our blockage analysis as compared to the reference 
model (1, Fig. 10). The relative amplitudes were computed and plotted in the same way as for Fig. 12. Lg blockage is also in this case moderate, with a 
maximum of about 9 dB. The additional Lg blockage of the perturbed Model 5 is generally less than 3 dB higher than for the corresponding 
unperturbed Model 2. 
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hence the current gap between observed and calculated Lg- 
blockage may in fact be less than presumed. The three sets of 
event records displayed in Figs 2-4 imply an Lg blockage of 
roughly 20, 0 and 6 dB, respectively, based on Ml magnitude 
differences using Pn and Lg amplitudes (Mendi & Husebye 
1995). By far the strongest blockage is obtained for Event 1, 
which also has the longest wave paths. For Event 2, there is no 
apparent Lg blockage for the coastal stations, except for MOL 
and NSS where the blockage is about 6 dB and takes place 
along the non-marine path segment. For Event 3, the Lg 
blockage is 6 dB on average. 

The North Sea grabens are obviously related to Lg blockage 
but not in a uniform manner. Our Event 2 (Fig. 3) is one such 
example, and Gregersen (1984) gives observational evidence 
of non-blocking Lg paths through North Sea graben areas. The 
Himalayas are also well known for efficient Lg blockage 
(Ruzaikin et al 1977), although non-blocking Lg paths can be 
found (Mykkeltveit & Husebye 1981). In other words, very 
efficient Lg blockage appears to be related to special crustal 
features like graben and mountain ranges. Strong intrinsic 
attenuation (low Q values) is one obvious candidate here while 
others may be sharp basin wedges towards crystalline rocks 
and extensive cracks and faults due to relative recent tectonic 
reworkings. Naturally, Lg-blockage mechanisms as mentioned 
here may also exist outside graben areas as exemplified by the 
Event 2 recordings at MOL and NSS. In a recent study, Xie & 
Lay (1993) demonstrated the problem of identifying clear Lg 
blocking mechanisms—not unexpected in view of the large 
fluctuations in the observational data. 

Finally, a possible limitation in our 2-D FD modelling 
experiments, besides that of not incorporating intrinsic Q 
attenuation, is that the code is for two dimensions. We 
intuitively expect to have more severe scattering in a 3-D 
environment and less scattered wavelets would return into the 
model ‘box\ On the other hand, grabens are dominantly 2-D 
tectonic features and most of the North Sea observational data 
are for paths nearly perpendicular to the elongated graben. 
There also seems to be a clear relationship between strong Sn 
and weak Lg, i.e. a considerable fraction of Lg energy leaking 
into the mantle reappears in the Sn coda. These features are not 
clear in our synthetics, as they would require a strong positive 
velocity gradient beneath the Moho and/or a more laminated 
or heterogeneous sub-Moho structure for proper modelling 
(Hestholm er a/. 1994). 


CONCLUSIONS 

In this study, we have simulated Lg blockage phenomena 
for realistic North Sea graben structures using 2-D finite- 
difference (FD) synthetics although intrinsic Q attenuation is 
not incorporated. The main results are as follows. 

(1) For extremely thin graben models with crystalline 
crust thicknesses of about 5 km plus sedimentary basins, Lg 
blockage amounts to 6-9 dB. 

(2) When an 8 per cent rms velocity perturbation was added 
in the graben crust only a few extra decibels Lg decay were 

(3) The dominant structural Lg-blockage mechanisms are 
the delay of the Lg waves in the thick sediments, Lg-to-Rg 
conversions (scattering) by lateral heterogeneities in the 
sediments, and Lg-wave leakage out of the crustal waveguide. 


(4) The above-mentioned Lg blockage mechanisms are 
more efficient in a thin crust than in a thick crust. 

(5) Part of the mantle waves are returned to the crust and 
appear in the Sn coda. 

(6) Real Lg observations require a somewhat strong block¬ 
age of the order of 15-20 dB which in turn require relatively 
strong intrinsic attenuation. A Q value of 100 at 2 Hz would 
suffice in this regard. 

(7) Our synthetics agree well with similar studies, although 
some authors concluded optimistically that synthetic blockage 
simulated real conditions well at ca. 6-10 dB amplitude 
reductions. 

(8) We also made comparisons with real observations 
(records from West Norway stations shown for three earth¬ 
quakes in the North Sea), where a strong Sn phase with an 
extended coda is often synonymous with strong Lg blockage. 

(10) Lg blockage has also been observed outside the North 
Sea Graben area, so blockage mechanisms must be complex and 
again low Q (strong attenuation) is likely to be a decisive factor. 

Our final concluding statement is that Lg exhibits very 
robust propagation characteristics and except for high intrinsic 
attenuation or very thin crustal waveguides in combinatin with 
thick sediments it is difficult to pinpoint specific blockage 
mechanisms. 
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Introduction 

The problem of artificial reflections from grid boundaries in the numerical discretization 
of elastic and acoustic wave equations has since long plagued geophysicists. Even if mod¬ 
em computers have made it possible to extend the synthetics over more wavelengths, 
equivalent to larger propagation distances, efficient absorption methods are still needed in 
order to minimize interference from unwanted reflections from the numerical grid bound¬ 
aries. In this study we examine applicabilities and stabilities of the Optimal Absorbing 
Boundary Condition (OABC) of Peng and Toksoz (1994, 1995) for 2-D and 3-D acoustic 
and elastic wave modeling. As a basis for comparison we use the Exponential Damping 
(ED) (Cerjan et. al. 1985), in which velocities and stresses are multiplied by progressively 
decreasing terms when approaching the boundaries of the numerical grid. 

Peng and Toksoz (1994, 1995) emphasized the importance of stability in the choice of 
aborbing boundary condition. As an example, they mentioned that Emerman and Stephen 
(1983) demonstrated that the boundary condition of Clayton and Engquist (1977) was 
unstable for a wide range of elastic parameters, and further that Mahrer (1986) found the 
boundary condition of Reynolds (1978) to be unstable as well. Peng and Toksoz (1994) 
gave some examples of 3-D elastic finite-difference simulations where OABC was used. 
The employed scheme was 4th order accurate in space and 2nd order accurate in time. 
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Elastic wave propagation in an unbounded homogeneous medium (without a free surface), 
were simulated, and the results were stable. They also made a comparison between the 
OABC, Reynolds’ (1978) and Higdon’s (1990) boundary conditions and found that using 
the OABC led to less artificial reflections than the two other schemes. Our study may be 
considered a further test of absorbing boundary conditions. We therefore implement both 
the OABC and ED in our 2-D finite-difference elastic wave algorithm (Hestholm and 
Ruud, 1994) and check their respective performances. Our goal is to find higher order 
accuracy algorithms using both methods. 

2-D finite-difference implementation 

We employ a 2-D finite-difference velocity-stress formulation (Hestholm and Ruud, 
1994), that solves the equations governing wave propagation in an elastic isotropic 
medium. Following Levander (1988) and Virieux (1986) we discretize the elastodynamic 
equations with two staggered numerical space differentiators. Details of this numerical 
discretization can be found in Kindelan et al. (1990), who developed optimal spatial finite 
difference methods based on the work of Holberg (1987). In the present work we use a 
method which is spatially accurate to 8th order in the interior of the computational 
domain. For time stepping a leap frog technique (accurate to 2nd order) is used. This 
allows us to achieve an upper bound of 1.5% for the relative error of the numerical group 
velocity at only 3 nodes per wavelength. A schematic of our staggered grid is given in Fig¬ 
ure 1. 

The OABC method extrapolates values on the numerical edges of a finite-difference 
grid. We express these values as a linear combination of the wave field at previous time 
steps and/or interior grids by exploiting the zeros and poles of the reflection coefficients in 
the complex plane. Approaching the grid edges, including the free surface, we apply suc¬ 
cessively lower order central, staggered finite-difference operators used for discretizing 
the spatial derivatives. We apply stress-free boundary conditions at the upper boundary. 
OABC is used at the bottom and sides of the grids. It is worthwhile to note that in the 
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implementation of the OABC it is important to avoid using comer points of the staggered 
grid. 

The ED method involves some restrictions on the useable part of our computational grid, 
since it is necessary to generate finite thickness damping strips along the grid edges. We 
exponentially damp velocity and stress in the strip by multiplying them by 

exp(-(a X A)^) . A is the number of grid points to the inner boundary of the strip, and 
a is a constant equal to 0.015. Damping is enforced within absorption strips of 20 grid 
points or 7 wavelengths at 60 Hz and 1 wavelength at 10 Hz along the bottom and the 
sides of the grid. For both OABC and ED, the coefficients are precomputed before time 
extrapolation is started. 


Test cases 

In our numerical tests the model size is 100 km by 70 km, or 40 by 30 wavelengths at the 
central frequency of 2.5 Hz, giving 351 grid points vertically and 501 grid points horizon¬ 
tally for a grid size of Ax = Az = 0.2 km. We design two different types of models in 
order to compare the artificial reflections using the OABC and the ED. The first model is 
homogeneous with a constant P-velocity of 7.1 km/s. The second model is a multilayer 
realization of the cmst and upper mantle consisting of several constant velocity layers 
(Figure 6). In all cases the S-wave layer velocity is related to the P-wave velocity via the 
Poisson ratio of v^/v^ = ^3. Densities p are linear functions of the P-velocities: 
p = 0.613 + 0.000328 • v^, given in kgim with given in m/s. The source is a Ricker 
wavelet (derivative of a Gaussian) in both space and time and is located 50 km from the 
left edge of the grid and 0.5 km below the surface. 

Stability of OABC 

We found that the OABC became unstable for long time simulations when using our 8th 
order accurate finite difference algorithm. We tried to use 8th order of accuracy in the inte¬ 
rior of the domain combined with second order along three levels adjacent to the absorb¬ 
ing boundaries (thus a very sharp discontinuity in the order of accuracy). We encountered 
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instabilities after 0(100) time steps, with the instability being initiated either along the 
absorbing boundaries or at the bottom comers of the computational grid. As already men¬ 
tioned, Peng and Toksoz (1994) achieved stable results for the OABC when used with spa¬ 
tially 4th order accurate schemes in a 3-D homogeneous medium. Thus for a realistic test 
of OABC using our scheme it was natural to experiment with lowering the spatial finite- 
difference order. This led to successful mns whether using 2nd or 4th order accurate finite- 
difference operators throughout the grid; in both cases we obtained completely stable mns 
with greater numerical dispersion. 

We therefore modified the OABC algorithm to achieve stability in our higher order 
scheme, thereby obtaining higher accuracy. We solved the problem by decreasing the 
finite-difference order gradually when approaching the boundaries. We had to perform 
several tests in order to find the proper number of grid points near the boundaries on which 
to apply a certain finite-difference order. It seems important to use 2nd order finite differ¬ 
ences at exactly two grid points adjacent to the boundaries. Inside this layer the number of 
grid points of higher order finite differences should be increased gradually with order up 
to where the 8th order method starts. 

In our implementation we arrived at a total number of 12 grid points adjacent to the 
boundaries at which the finite-difference order were gradually lowered. We used 8th order 
in the interior of the grid, then 6th, 4th and 2nd order as approaching the boundaries. The 
following procedure was found to be optimal (in the sense that it delayed instabilities the 
longest) in our scheme when approaching the grid boundaries: 6 layers of 6th order finite 
differences, then 4 layers of 4th order and finally 2 layers of 2nd order adjacent to the grid 
boundary. In this way we were able to ran for at least 14000 time steps (about 60 seconds 
of simulated time) before a slowly growing instability started appearing along the absorb¬ 
ing boundaries. Simulating wavefields for such large time laps enable the most interesting 
wave phenornena to be adequately synthesized. 

The OABC requires the zeros and poles of the reflection coefficients to be in the first 
quadrant of the complex plane (Peng and Toksoz, 1994,1995), which leads to < 7t/2, 
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where A;c is the grid size and k the wavenumber. From this it follows that Ax < %/ 4, 
where X = 1% Ik is the shortest wavelength. Therefore the OABC requires at least 4 sam¬ 
ples per shortest wavelength to ensure stability, at least near the grid edges. In our context 
of a broadband Gaussian point source this cannot be ensured, since we do not have explicit 
control over the complete frequency range that propagates. The Ricker wavelet source is 
widely used in seismic simulations. However, the central frequency of the source is 2.5 
Hz, which leads to a wavelength of 2.8 km for our choice of P-wave velocity in the homo¬ 
geneous medium. This leads to 14.2 samples per wavelength, which is 3.5 times more than 
is needed for the OABC. Nevertheless, the source excites also a high-frequency part that 
seems to destroy stability of the OABC when used together with higher order finite differ¬ 
ences. The high-frequency part generally needs finer sampling in order for the stability cri¬ 
terion to be satisfied, and instability occurs specifically in connection with higher order 
finite differences. Of course, an obvious solution to this problem would be to replace the 
Ricker wavelet with a bandlimited wavelet. 

Another aspect that should be emphasized is the use of a free surface in our applications. 
This leads to Rayleigh waves propagating with much smaller speeds and wavelengths than 
the P- and S-waves, which consequently may violate the stability criterion of the OABC. 
This problem is common with absorbing boundary conditions, and the effect becomes 
even more prominent in the case of free surface topography. 

For reasons described, and particularly in order to be able to use realistic sources in con¬ 
nection with OABC, we developed the given optimal procedure of implementing OABC 
in our higher order finite-difference scheme. Following this procedure, instability is 
delayed until the most interesting wave phenomena have been synthesized. 

Efficiency 

In the first simulations the homogeneous medium model was used, using (1) the OABC 
(Figures 2 and 3) and (2) the ED (Figures 4 and 5). All plots were scaled according to the 
maximum value anywhere in the grid at each time step. If a permanent maximum value 
were used, reflections would be invisible. In Figure 2 the -norm of the velocity vector is 
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displayed in a snapshot 11s after a pressure wave has been initiated. Here, the P-wave 
front has just reached the bottom of the model. The snapshot reveals a weak P-P reflection 
from the bottom, while stronger artificial P-P and P-S reflections can be seen from both 
sidewalls. Interestingly, the reflections from the left sidewall are stronger than the ones 
from the right sidewall. The reason is our staggered grid description of Figure 1. Since we 
use the same number of grid points of each field variable in each grid dimension, a differ¬ 
ent set of variables generates the left and right grid edges (Figure 1). Specifically, the three 
(of totally five) variables 8^^, 8^^ and w are absorbed further inside the grid on the right 
than on the left grid edge and explains the better absorption on the right edge (Figure 2). 
Figure 2 also shows a P-SV head-wave propagating from the upper left comer only. This is 
because the velocity analog of the stress-free surface conditions (Hestholm and Ruud, 
1994) is implemented in an asymmetric way according to Figure 1. 

In Figure 3 the velocity vector modulus (L 2 -norm) is displayed in a snapshot 22 s after 
the pressure wave has been initiated. Predominantly artificial reflections can be identified 
here, since the dominant non-artificial wavefield has passed out of the model frame. The 
reflections from the S-wave are stronger than from the P-wave due to our choice of maxi¬ 
mum P-wave absorption (Peng and Toksoz, 1994,1995). The strongest reflections are seen 
to come from the bottom, while the vertically oriented wave front seen near the right 
boundary is the P-wave reflection from the left boundary. 

Figures 4 and 5 are similar to Figures 2 and 3, except that ED was applied within 20 grid 
points wide strips along the numerical boundaries at the bottom and sides. Clearly, the ED 
scheme works better than the OABC in absorbing both P-waves and S-waves. It absorbs 
equally well along the bottom and the sidewalls. In Figure 5 the non-artificial waves have 
passed out of the model, and only artificial grid reflections remain. Investigating amplitude 
values at the layer adjacent to the grid edges, these reflections were found to be weaker 
than the corresponding reflections with the OABC by a factor of approximately 10“^ , and 
furthermore, the run is completely stable. Again, because of the staggering, a Rayleigh 
wave can be distinguished propagating along the surface on the right part of Figure 5, 
whereas the one on the left is too weak to be seen. 
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Application of lower order finite differences near boundaries introduces extra disper¬ 
sion. This dispersion seems to be less well handled by the OABC than by the ED. The rea¬ 
son is probably that high frequencies are quickly attenuated by ED. Methods like ED 
(Cerjan et. al., 1985) are known to absorb reflections of relatively short wavelengths better 
than most methods. Because of this, we see that ED absorbs both the high-frequency part 
from the broadband source and the Rayleigh waves at the surface quite well. In our exam¬ 
ples, the high-frequency part from the source seems especially problematic for the OABC 
to absorb. The problem may be partly due to unspecified behavior of the OABC when it is 
used to absorb frequency parts for which its stability criterion is not satisfied. 

A homogeneous model is an oversimplification of the real Earth, and so further compar¬ 
ison between OABC and ED is tied to wavefield simulation using a multilayered crust - 
upper mantle model shown in Figure 6. Velocity vector modulii (amplitudes) as in Figures 
2-5 are shown in Figures 7 and 8, but now using the background model of Figure 6. This 
fact, together with the multiple additional reflections from the various layer boundaries, 
cause these snapshots to be less clear than those for the homogeneous medium model 
(Figures 2-5). Figure 7 shows the wavefield 22 seconds after the pressure wave was initi¬ 
ated, using OABC. Also here asymmetric P-P and P-S reflections from the sidewalls can 
be seen, together with long term effects of wave conversions and reflections from the lay¬ 
ers. The absorption at the bottom is seen to be somewhat better than that at the sidewalls 
due to the staggered grid. 

Figure 8 shows the corresponding run using ED. The absorption strips of 20 grid points’ 
width are clearly visible along the numerical boundaries. In comparing amplitudes at loca¬ 
tions of artificial reflections from the last two runs we found them on the average 0(10) 
stronger from using OABC than from using ED. 

Conclusions 

While a naive implementation of the OABC in our 2-D finite-difference scheme resulted 
in instabilities, a simple modification of the original algorithm allowed us to delay them 
beyond the time of interest. We arrived at an optimal procedure for implementing Peng 
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and Toksoz’s OABC algorithm, and then compared the OABC results with those obtained 
by using ED. 

In implementing the OABC we found the choice of spatial finite-difference order to be 
important. With broadband Gaussian point source types the stability criterion of the 
OABC is not fulfilled for the complete frequency range, even though stable behavior is 
achieved by using 2nd and 4th order spatial finite differences throughout the grid. OABC 
always behaves unstably when using higher order spatial finite differences. The order of 
accuracy must be decreased when approaching the numerical boundaries in order to delay 
the instability. In our synthetic discretization scheme an 8th order accurate spatial finite- 
difference approximation is used in the interior of the grid. In moving from the interior 
towards the boundaries we found the following procedure to be optimal: 6 layers of 6th 
order finite differences, then 4 layers of 4th order and finally 2 layers of 2nd order adjacent 
to the numerical boundary. It seems clear that ED absorbs better than OABC. In a homo¬ 
geneous medium, the artificial reflection amplitudes from using ED within a thin strip 
were smaller by a factor 10"^ compared to using OABC. A reason might be insufficient 
sampling of higher frequencies near boundaries using OABC, and thereby violating its 
stability criterion for this frequency range. Another reason is the OABC’s apparently 
smaller ability to handle dispersions arising from using lower order finite differences near 
boundaries. The mentioned drawback of additional computer storage are offset by the 
shorter simulation time using ED. In our 2-D applications we absorb at 13.2 % of the total 
number of grid points. In a 3-D application using the same grid as here with the second 
horizontal dimension length equal to the first, we would use 6.3 % of the grid points for 
absorption. 
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Figure Captions 

FIG. 1. Schematic of staggered grid variable definition. 

FIG. 2. Snapshot of the velocity vector amplitude at t=l 1 s after a pressure wave has been 
initiated from a Gaussian point source at depth 0.5 km below the surface. The medium is 
homogeneous and OABC is used. 

FIG. 3. Snapshot of the velocity vector amplitude for the same situation as in Figure 2 at 
t=22 s. 

FIG. 4. Snapshot of the velocity vector amplitude for the same situation as in Figure 2, but 
in this case ED is used. 

FIG. 5. Snapshot of the velocity vector amplitude for the same situation as in Figure 4 at 
t=22s. 

FIG. 6. Snapshot of the velocity vector amplitude at t=1.2 s after a pressure wave has been 
initiated from a Gaussian point source at depth 0.5 km below the surface. The source posi¬ 
tion and the background model of the layered medium are clearly visible. The maximum 
P-velocity is 8.18 km/s right below the Moho. The Moho is located at 35 km depth. The 
closest layer above the Moho represents a P-velocity of 6.9 km/s. The next layer is the 
maximum velocity layer in the crust with a P-velocity of 7.1 km/s. Just below the free 
plane surface there is the minimum velocity layer with a P-velocity of 6.6 km/s. 

FIG. 7. Snapshot of the velocity vector amplitude at t=22 s after a pressure wave has been 
initiated from a Gaussian point source at depth 0.5 km below the surface. The medium 
model is layered (Figure 6) and OABC is used. 

FIG. 8. Snapshot of the velocity vector amplitude for the same situation as in Figure 7, but 
in this case ED is used. 
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ABSTRACT 

Three-dimensional (3-D) finite-difference (F-D) modeling of seismic scattering 
from free surface topography has been pursued. We have developed exact 3-D free 
surface topography boundary conditions for the particle velocities. A velocity-stress 
formulation of the full elastic wave equations together with the boundary conditions 
has been numerically modelled by an 8th order F-D method on a staggered grid. 
We give a numerical stability criterion for combining the boundary conditions with 
curved grid wave equations, where a curved grid represents the physical medium 
with topography. Implementation of this stability criterion stops instabilities from 
arising in areas of steep and rough topographies. We have simulated scattering from 
teleseismic P-waves using a plane, vertically incident wavefront and real topography 
from a 40 X 40 km area centered at the NO RES S array of seismic receiver stations 
in south-eastern Norway. Synthetic snapshots and seismograms of the wavefield show 
clear conversion from P— to Rg— (short period fundamental mode Rayleigh) waves in 
an area of rough topography approximately 10 km east of NORESS. This result is 
consistent with numerous observations. 


1 



INTRODUCTION 


The theory presented here is a direct extension from 2-D to 3-D of the corre¬ 
sponding theory in Hestholm and Ruud (1994). Inclusion of topography at the free 
surface of an elastic medium leads to improved modeling of near-surface scattering 
effects, especially those in the high frequency part of the wave field. Relatively little 
has been published on the modeling of free surface topography in 2-D and even less on 
the modeling of 3-D surface topography. This applies both to F-D methods and any 
other numerical discretization method. However, a work on this problem is that of 
Frankel and Leith (1992) who used an F-D scheme of fourth order accuracy in space 
and a density taper to zero starting at the height of the free surface while keeping 
the medium P—velocity unaltered. They achieved stable results by multiplying the 
crustal density by 0.4 for the locations one grid point above the free surface, by 0.16 
for the locations two grid points above the free surface, and so on. In this manner 
they were able to obtain reasonable modeling results for free 3-D surface topography 
models. 

More recently, Tessmer and Kosloff expanded their procedure for elastic wave 
modeling with free surface topography from 2-D to 3-D (Tessmer et al., 1992; Tessmer 
and Kosloff, 1994). They use a spectral discretization in space, being different in the 
horizontal and vertical directions. The velocity—stress formulation of the elastic wave 
equations is transformed into a curved grid. At the free surface, the stresses and 
velocities are transformed into a local system in which the vertical coordinate axis 
is normal to the surface element. The free surface conditions are then implemented 
by a ’characteristic’ treatment of the velocity and stress components, before rotated 
back to the original system. In this method, both velocities and stresses are rotated 
into the local system at each point of the topography surface. 

In the present work explicit 3-D boundary conditions for a free surface topogra¬ 
phy are derived. The basis is the vanishing stress condition for a free surface. As 
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in Tessmer and Kosloff (1994), a 3-D curved grid which is stretched in the vertical 
direction is adapted to the surface topography, i. e. the top surface of the grid co¬ 
incides with the surface topography. A coordinate transform is then introduced for 
transferring the elastic, isotropic wave equations from the curved to a rectangular 
grid in which the numerical computations are done. At the topography surface, the 
velocity boundary conditions for a free surface are implemented into a local, rotated 
system at each point on the surface. Each of these systems has its vertical coordi¬ 
nate direction coinciding with the direction of the normal vector of the surface at 
the given point. The velocity boundary conditions are subsequently rotated back to 
the rectangular system. Once the boundary conditions are given in this system, the 
n um erical discretization can be performed. 

In the following paragraphs, the 3-D equivalents to the 2-D equations (Hestholm 
and Ruud, 1994) will be stated. In each instance, it is possible to verify coincidence 
between the 3-D and the 2-D case by eliminating the r-direction (2nd coordinate 
direction of the computational grid) in 3-D. We can also verify the coincidence be¬ 
tween the plane surface conditions in 3-D and the 3-D surface topography conditions 
applied to a plane surface. A description of the numerical discretization will be given 
and stability criteria for the method will be assessed. Then we present simulated scat¬ 
tering from teleseismic P—waves using a plane vertically incident wavefront and real 
topography from an area centered at the NORESS array in south-eastern Norway. 
Snapshots and synthetic seismograms of the wavefield will be shown, clearly display¬ 
ing Rg waves in areas of rough topography. Finally, we look at future prospects for 
3-D F-D modeling, particularly in light of the recent parallellization of the seismic 
code. 


3 



ELASTIC WAVE MODELING FORMULATION 


The basic equations governing wave propagation in a continuous elastic medium 
are the equations of motion and the stress—strain relationship. The velocity—stress 
formulation (Achenbach, 1975; Virieux, 1986) can be written in 3-D as 


XX , dcfxy , 9^xz j f 
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OX ay az 

d<Txy j d^yy , dcTyz ^ ^ 
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where p is the density and A and p axe Lame’s parameters, fx, fy and fz are the 
components of the body forces, u, v and w are the particle velocity components and 
(Txx, CTyy, (Tzz, CTxy, (^xz and (Tyz are the stress components. 

We introduce a linear mapping from a rectangular (^, r, 7?)-system (Figure 1) 
to a curved {x,y,z)-system (Figure 2), where both systems have positive direction 
upwards for the vertical coordinate. The 3—D mapping can be written as 
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where 2 ;o(^,t) is the topography function, and the rectangular (^, r, 7j)-system is 
limited by ^ = 0, ^ = ^max, r = 0, r = Tmax, tj = 0 and rj = r/max- For the 
curved {x,y,z)-system the degree of stretching is proportional to the distance from 
the bottom plane of the system {z = 0). From equations (10)-(12) we get, for a 
differentiable function f{x, y, z), 

^ = (13) 

dx dy dx ’ 

+ (14) 

dy dr dydy' 

dz dydz' ^ ^ 

Expressions for the partial derivatives, which are needed in the medium equations, 
are found from equations (10)-(12) (see Appendix A), 

g-. t- 
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The velocity-stress formulation of the equations of motion and Hooke’s law is 
given in the curved (x, y, z)-grid by equations (l)-(9). Expanding these by the chain 
rule (Appendix B) as in Hestholm and Ruud (1994), we get the medium equations in 
the rectangular (^, r, ? 7 )“System, 
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Equations (21)-(29) are the momentum conservation equations and Hooke’s law given 
in the rectangular (if, t, T/)-system. 


FREE SURFACE BOUNDARY CONDITIONS 


The 3-D free boundary conditions for the velocities at a locally horizontal surface 
(or in a system where the z—axis is parallell to the local normal vector of the surface) 
resulting from the vanishing stress condition can be written 


du _ dw 

dz dx ’ 
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dz dy ’ 
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(32) 


with X and y the horizontal coordinates and 2 the vertical coordinate. We want to 
apply these conditions to a topography surface. At each surface point, we introduce a 
local coordinate system [x', y', z') in which the z'-axis coincides with the local normal 
vector direction of the surface. In this local system we impose the conditions (30)- 
(32). Once this is done, they have to be rotated back to the rectangular system 
((^, T, rj) where the numerical computations can be done. This rotation is expressed 

by 


V = Av 


(33) 
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where v and v ' are the particle velocity vectors in the (^, r, r\)- and the (x', y', z')- 
systems respectively. A is the rotation matrix, which can be given by 

/ cos 0 0 sin0 \ 

A = — sin 0 sin ^ cos (j) cos ^ sin • (34) 

^ — sin 6 cos <f) — sin (f) cos 6 cos (p / 

where 0 is the rotation angle between the ^-axis and the local x'-axis in the (^, ??)- 
plane and is the rotation angle between the r-axis and the local y -axis in the 
local (y', 2 :')-plane. We also have the relations tan0 = dzo{^,T)fdC and tan<ji = 
dzo{^,T)/dT COS0 (see Appendix C). 

Now the calculations of the rotation from the local (x', y', . 2 r')-system back to the 
(if j r, y)—system can be performed as in Appendix C. We arrive at the 3—D boundary 
conditions for free surface topography given in the computational (^, r, y)-system by 
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e = cos [arctan (d)], 

(40) 

/ = sin [arctan (d)], 

(41) 

ll 

(42) 

q = cos [arctan (p)], 

(43) 

r — sin [arctan (p)]. 

(44) 


Equations (35)-(44) are exact 3-D boundary conditions for an arbitrary, smooth, 
free surface topography. They result from rotating the velocity free surface conditions 
from local systems at each point of the surface topography into a rectangular system 
(Appendix C). The vertical axis in each of the local systems is normal to the local 
topography. Also, the boundary conditions (35)—(44) are obviously not restricted to 
the F-D method or any other numerical discretization technique. 

NUMERICAL DISCRETIZATION 

For the numerical discretization, we refer to the corresponding paragraph in Hes- 
tholm and Ruud (1994). The same spatial and time discretization methods as in 
the 2-D case (Kindelan et ah, 1990) are used. The schemes employ a staggered 
discretization stencil as in Levander (1988) of the velocity-stress formulation of the 
elastodynamic wave equations, which is again based on the work of Virieux (1986). 
An advantage of using a staggered definition of variables is that we can avoid ex¬ 
plicit definition of the stresses at the surface topography as it suffices to define the 
velocities there. In order to get the velocities and stresses explicitly defined at each 
time step, we have to stagger the vertical velocity component w one half grid length 
downwards. Corresponding numerical definitions of the other variables have to be im¬ 
plemented. Generally, u is staggered one half grid length in the positive ^-direction, 
V is staggered one half grid length in the positive r-direction, while w is staggered 
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one half grid length in the negative r?-direction (downwards). The 3-D boundary 
conditions, equations (35)-(37) are discretized by 2nd order, staggered F-D opera¬ 
tors (Fornberg, 1988). Below the free surface, the central, staggered F-D method’s 
order is gradually increased with depth, via 4th and 6th up to 8th order, the latter 
is the order used inside the domain. This 8th order method is dispersion-bounded 
and cost-optimized (Kindelan et ah, 1990). Along the artificial grid boundaries ex¬ 
ponential damping according to Cerjan et al. (1985) is used. In a layer of 20 points 
along each grid boundary, the stresses and velocities are multiplied by exponentially 
decreasing terms towards the boundary. 

To find the velocity components at the surface topography from the closed system 
(35)-(37), we solve it directly as a linear system with respect to the velocities at the 
surface as they are defined in the vertical derivative discretizations. In this procedure, 
the horizontal derivatives are defined half a grid length below the free surface, and 
the derivatives are discretized by 2nd order finite differences. 

STABILITY CRITERION FOR THE SURFACE TOPOGRAPHY 

MODELING 

A numerical stability criterion for the surface topography modeling comes from 
the equations of motion in the medium. A necessary condition to keep a run stable 
is that the absolute values of the parameters A{^,T,ri), B{(,T,ri) and C{^,t) all be 
kept less than 1. This necessary condition can be written as 

' 1 , 

> max ^ ||^^o(6r)/^^||, ^ 

'^max 

. ||52ro(^,T)/5r||, 

where r]max is the total depth of the numerical (^, r, 77)-grid and 52o(^,r)/5^ and 
dzo{^-,T)ldT are topography slopes. This condition must be satisfied at every point 
on the surface. Note that this condition can always be satisfied, if necessary, by 
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uniformly increasing everywhere. Equivalently, rj^nax can be decreased to 

satisfy the inequality. 

For some combinations of sources and reliefs, this condition might not be suffi¬ 
cient to avoid an unstable growth right below the surface at the steepest parts of 
the topography. In these cases, rjmax should be reduced enough for stability to be 
maintained. By decreasing r^jriax rather than increasing uniformly, the orig¬ 

inal physical model (represented by the curved grid) is conserved. It is a matter of 
experimentation, then, how much rjmax niust be reduced in order to make wavefield 
simulations stable for a specific relief. The number of vertical samples must be the 
same in the curved and rectangular systems. This means that in order to maintain 
accuracy and simultaneously retain stability gained from increased curved to rect¬ 
angular grid depth ratio, the vertical grid distance (vertical distance between grid 
points) in the numerical (^, r, ? 7 )-system must be reduced by an appropriate factor. 

For most sources/topography constellations, it turns out that min{{zo{^jT)} ^ 
^max is enough for a simulation to be stable. Therefore, this criterion might be stated 
as a sufficient condition for stability in most cases. From our experiments, this order 
of magnitude for the ratio between the physical model and the numerical grid depth 
appears necessary whenever the topography data exhibits rough behaviour (large 
spatial second derivatives) near its steepest slopes. For rough topography without 
large spatial derivatives or steep topography without large second derivatives, this 
ratio may be smaller in order to achieve stability. In many such cases the stated 
necessary stability condition is also sufficient. 

P~ TO RG-SCATTERING FROM TOPOGRAPHIC RELIEF 

Since 3-D topography surfaces in general, and real reliefs in particular, are sel¬ 
dom used in detailed wavefield modeling, we hope that the above approach should 
constitute a powerful method for the F-D elastic wave modeling including surface 
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topography. The main goal of the research performed is to improve our understand¬ 
ing of wave propagation and scattering phenomena for local and regional distance 
travel paths. An interesting case here is observationally well documented. Scatter¬ 
ing of teleseismic P-wave energy into Rg have been extensively studied using data 
from the 3 km aperture NORESS array in south-eastern Norway (Bannister et ah, 
1990; Gupta et ah, 1993; Hedlin et ah, 1991; Hedlin et ah, 1994). For our 3-D F-D 
simulation of this phenomenon we have obtained digital elevation data for an area 
of 100 X 100 km centered on the NORESS array (Figure 3). Due to the huge com¬ 
puter memory requirements of 3—D F—D methods, we have so far been restricted to 
a model of size 40 x 40 x 35 km with 0.2 km sampling. However, this situation is 
now improved by parallellization of our code. In all the examples shown the incom¬ 
ing wave is a vertically incident plane P—wave simulating a teleseismic short period 
P~phase. The center frequency of the Ricker wavelet is 2.5 Hz, the P-wave velocity 
of the homogeneous medium is 6.0 km/s and the Poisson ratio is 0.25. 

In our present experiment the elevations were multiplied by 0.5 so that topography 
undulations were damped. By doing this the magnitudes of A((^, r, 77 ) and 5(^, r, 77 ) 
were reduced so as to achieve stable results without the need of increasing the model 
depth beyond the depth of the numerical grid. Alternatively, the numerical grid depth 
could have been reduced or the model grid depth could have been increased by some 
factor. Snapshots of the wavefield (vertical component of the particle velocity) are 
shown in Figure 4. A dominant feature here is the low-frequency/long-wavelength 
artificial reflections from the absorbing model boundaries. Although the exponential 
damping technique is the most efficient absorbing boundary method we have tested, 
a wave incident at an angle of 90® with the boundary is the most difficult case. 
Fortunately, the artificial reflections are frequency dependent, with longer wavelengths 
than most of the surface waves scattered by the topography, and can therefore be 
removed by filtering and image processing techniques. As seen from Figure 4 the 
scattered surface waves appear to radiate out from secondary point sources which 
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coincide with areas of high topographic gradients (Figure 3). The dominant scattering 
points are along the steep valley side east of ’Bronkeberget’ about 10 km east of 
NORESS (Bannister et ah, 1990). Also in seismograms from the west-east profile 
P-to-Rg scattering from this area is clearly seen (Figure 5). 

In order to quantify and locate areas of significant topographic features, we have 
computed different functions depending on first and second order derivatives of the 
topography function z — ZQ{x^y). First, we computed the slope, i.e., the length of 
the topography gradient vector as 


|V2o(x,y)| 



(45) 


This vector is useful for defining areas of strong scattering. Additionally, smooth 
topography is assumed in deduction of the boundary conditions. As a measure for 
topography roughness, we compute the Frobenius norm of the Hessian matrix of 
zo{x,y), i.e.. 


\\H{x,y)\\F^ 
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Both these functions have their largest values in the areas of strongest scattering 
(Figure 6). 

Recently, we modified the code to run on parallel computers. The tool used for 
parallellization was MPI (Message Passing Interface), which is the first attempt of 
standardization to some degree of parallel programming tools. Besides speeding up 
computations, the parallellization also allows us to use a larger computer memory, 
which is essential for realistic seismic wavefield simulation. Test runs employing up 


to 2GBytes of memory have been performed. For a simulation of this size we used 48 
hours on 48 processors on the Intel Paragon machine at the Univ. of Bergen, Norway. 
Increasing the number of processors to 80, which is the maximum possible number of 
processors on the machine’s batch queue, the simulation time was reduced to a little 
more than 6 hours, i. e. a speed-up factor of 8. The main reason for the significant 
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speed-up in the latter case is that the complete 2GB of memory was distributed over 
real memory on the processors, so that no virtual memory had to be used, and thus 
paging to and from disk was avoided. 

DISCUSSION AND CONCLUSIONS 

Even for small scale experiments 3—D F-D synthetic wavefield analysis provides, 
as demonstrated here, improved insight and a better understanding of surface scatter¬ 
ing phenomena. We find it particularly gratifying that the strong P-to-Rg scattering 
from the Bronkeberget hill, observed by Bannister et al. (1990) through analysis of 
NORESS recordings, can be realistically synthesized. Another interesting feature is 
that for certain spatial low-frequency variants of the local NORESS topography the 
Rg-wave propagation seems to change abruptly over relatively small distances. Prob¬ 
ably, this is due to strong directionality dependence of the scattering from some topo¬ 
graphic features. Such wavefield characteristics are sometimes observed in NORESS 
recordings and at the German array GERESS located in the Bavarian hills - at some 
sensors the Rg-phase is prominent while hardly visible at nearby sensors less than a 
kilometer away. In other words, also small scale crustal features may contribute to 
blocking effects as often observed for Rg- and Lg-phases across structural obstacles. 
Many of such wavefield phenomena are out of range for 2-D F-D synthetic experi¬ 
ments; hence our emphasis on continuous research efforts on 3—D synthetic simulation 
of crustal wave field propagation. 

Future directions of our research will be to increase the model size to also include 
the scattering areas south-east of Lake Mjpsa (Bannister et ah, 1990). Furthermore, 
in order to compare directly with NORESS recordings it is necessary to allow for 
non-vertical incident waves and to use source time functions derived from teleseismic 
P-beams. We also plan to implement the topography boundary conditions into 2-D 
and 3-D viscoelastic schemes. This will be done in order to be able to simulate a 
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realistic attenuation of waves, in particular the strong waves generated from areas 
of significant surface topography. In addition, this will greatly improve the absorb¬ 
ing boundary conditions. Although the physical dimensions of 3—D models used in 
simulation will still be relatively small compared to 2-D models employed, the use 
of parallel computer technology opens new avenues to the study of 3-D scattering 
phenomena in fields like random media and corrugated interface scattering. 
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APPENDIX A-PARTIAL DERIVATIVES IN MEDIUM EQUATIONS 


For the medium equations we need drjldx, drfjdy and dyldz. They are found 
from 


dx dx dr dx dy 

d^ dx dr dx'^ drj dx ’ 

(A-1) 

_ Q 

d^ dx dr dx dy dx ’ 

(A-2) 

dz d^ dz dr dz dy 
d^ dx dr dx'^ dy dx ’ 
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This leads to 
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where 

dxfdydz dzdy\ dxfdydz dy dz \ dx f dydz dz dy \ 

d^\dTdri drdr]) dT\d^d'ri drjd^J dr]\d^dr d^dr) 

With our choice of mapping functions, equations (10)-(12), we get 
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From this we get the expressions (16)-(20). 


(A-17) 
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APPENDIX B-MEDIUM EQUATIONS 


Applying the chain rule to equations (l)-(9) and using the properties of equations 
(13)-(15) leads to 


du _ da^x , dcTxx dri da^y da^y dr] ^ XZ drj 

^ dt ~ dr] dx dr drj dy dy dz 

^ _ daxy daxy drj d^ dcTyy dy d^^ , 

^ dt d^ drj dx dr drj dy drj dz 

^ _ d^ d^^ d^ d^^ d^^ 

^ dt ~ d^ drj dx dr drj dy drj dz 

dcTxx _ y dudrj ^ I 2 (— + 

dt ~ \d('^ drj dx^ dr drj dy drj dz) ^ \ d^ drj dx) ’ 

da-yy _ s (9u dudrj ^ ^drj\ f — 4 . 

dt ~ \d^^drjdx dr drjdy drj dz) ^ \dT drjdy)' 

dcTzz \ f dv drj dy^dj^ / dw drj \ 

dt ~ \d^'^ drj dx^ dr drjdy drj dz) ^\drjdz)' 

doxy f dv dv drj du du 
dt ~ ^ \d^ drj dx ^ dr drjdy) ' 
duxz _ (dw dwdrj du 
dt ^ \d^ drj dx drj dz) ^ 
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(B-5) 

(B-6) 

(B-7) 

(B-8) 

(B-9) 


Substituting for drjjdx, drjfdy and drj/dz from equations (16)-(20), we get the 
medium equations (21)-(29). 
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APPENDIX C-BOUNDARY CONDITIONS 


Assume a velocity vector v with components u, v and w is given in the (^, t, r/)- 
coordinate system with basis vectors i, j and k. This system is then rotated through 
angles into a new (a:', 2 :^)-coordinate system with basis vectors i', and fc'. 

9 is the rotation angle between the ^-axis and the a;'-axis in the ((^, y)-plane. (j) is the 
rotation angle between the r-axis and the j/'-axis in the {y\ ^')-plane. In this new 
system the vector v is denoted by u' with components u\ u' and w'. Then we have 
the relationships 

/P\ /i\ 

I' = A j , (C-1) 

\k') \k) 

where A is the rotation matrix, given by equation (34). Correspondingly, 

Z \ ( \ ^ 

j =A-^ /' =A^ , (C-2) 

yk) \k') \k'J 

where A~^ and A^ are the inverse and transposed of A respectively. 

Using |n| = ]J + 1, a unit normal vector to a surface topography 
element can be written as 


- n 1 / dzo{^,T) dzo{^,T) 

^ “ Ini “ Ini I ’ ar ’ 


= (— cos (f) sin 0, — sin (f>, cos ^ cos 6 ) 


with our choice of rotation angles. From this we get 
tan<^ dzoU,T) . 

-^ =i.e. tan^ - '-cosO, 

cos 9 OT OT 

J \ ,. f r) ^ 

cos <p = cos arctan I — -cos u 

and sin (p = sin arctan f ——- cos 9 
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The coordinate transformation for v is given by v = Au', or v' = A~ v. Compo¬ 
nentwise this is 


u' = (cos 6 )u — (sin 6 sin 4 >)v — (sin 6 cos 4 >)w, 

v' — (cos (f>)v — (sin 4 >)w, 

w' = (sin 9)u -\- (cos 6 sin </>)?; + (cos 6 cos 

Applying the chain rule to a differentiable function /, we then get 

^ + 8/1!* = K cos « + sin 

dx' dx' dr dx' dr] dx' d^ dr] 

df ^ 

dy' d^ dy' dr dy' drj dy' 

d f ^ J" ^ f 

= sin 6 sin 6 ) + — cos <f> -I- cos 9 sin</>, 

d^ dr drj 

df ^ 

dz' d( dz' dr dz' dy dz' 

— ^(—sm9cos4>) + ^(—sm<f>) + cos 9 cos f. 
dC dr^ dy 


(C-IO) 


(C-11) 


(C-12) 


The free boundary conditions (30)—(32) for the velocities have to be enforced in the 
local {x',y',z')-system, where the ^-axis is normal to the surface at the local point, 


du' dw' 

dz' dx' ’ 

dv' dw' 

dz' dy' ’ 

dw' A fdu' dv' 


dz' \-\-2y,\dx' dy')' 

If the chain rule is applied as above, we get 

du'. du'. . . du' , 

-—(—sin^cos^) + -^(—smf) + — cost^cos^ 

d^ dr dy 

dw' dw' 

= — —cos9- —siii9, 
di dy 

dv', . . i\ , L a A 

sin9 cos <p) + -;;^(—sincp) + cos cos <p 

d^ dr dy 


(C-13) 

(C-14) 

(C-15) 


(C-16) 


dw' ■ n ■ J.\ A 0 ■ A (C 

- — 7 ^{- sm 9 sm A) - cos f —r— cos 9 sin (p, (C-17j 

d^ dr oy 


21 



dw' - dw' . ,x a A 

-—(—smOcosS) + -p—[—sm(p) + -r—cos cos <p 
OT or] 


A (du' du' . .,dv' . • ,x 

= -XTv (,W “ a? 


dv dv A ■ I \ 

+ -r— cos <P + - 7 ^ COS 8 sm<p]. 
OT Or] 


(C-18) 


Now we apply the above expressions for u', v' and w' that were obtained from the 


rotation. This leads to 


+ 

+ 


+ 


+ 

+ 


.du . . dv . dw\ a ,x 

cos^— — smpsmfflTT- — smpcos^Tr- ( — smv cos 9 ) 

di di d^ 

.du . . . ,dv . dw] . , 

cos 07 ^-sin^sincpTr-smycos<p^ (—sm^j 

Or dr Ot _ 

.du . n . idv . . ,dw] , 

cos0-^ -sin 0 sm * 7 ;-sin 0 cos <p-;— cos 0 cos (p 

dr] dr] Or] 

• Adu . . ,dv 1.9'^] A 

— sm0— — cos0sm<p— — cos0cosp-^ cos0 

d^ d^ d^_ 

■ a9u a ■ ,dv . J^dw] . 

— sm 6~ -cos u sin 6— -cos u cos cp-^ sin c/, 

dr] dr] dr] 

dv . dw] a A\ 

cos<;d^ — sinip-^ (—sint'cos^j 

dv . dw] , . ,x 

cosw-r -sinffiTT- (—sin 9 ) 

dr dr _ 

,dv . .dw] 

cos <p— -sm cos u cos <p 

or} OTj 


. du dv dw 

— sm 0— — cos 0 sm (p— — cos 0 cos p— 

d( di d^ 


+ 

+ 


+ 

+ 


• a ■ n 

— sm u-T -cos d sm o— -cos u cos 9 -^;— 

dr or or 

. du . . ,dv dw 

--sm0— -cos 0 sm p— -cos 0 cos p-^ 

dr] dr] dr] 

. Adu . . ,dv dw 

sm fc'Trr + cos 0 sm p— + cos 0 cos p— 

d( d( 

. „du . . ,dv _ . ,dw 

sm 6——h cos d sm (p-p :—h cos v cos 

or or or 


(— sin ^ sin (j)) 
cos (f) 

cos 6 sin </>, 

(—sin^cos (f) 
(— sin <f>) 


. ^ ^ ,dv . ..dw 

sm6— + cos 6 sm 0— + cos d cos 0—, 
dr) dr] orj J 


A -}- 2/i 


cos 6 cos 0 
dw 


.du . . dv . 

cos a— — sm d sm 0— — sm tl cos 0 


cos 9 


(C-19) 


(C-20) 
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+ 

+ 

+ 

+ 


du . dv . dw 

cos 0-r -Sin sm 4>— -sin d cos <p-^ 

or] 


dv . dw 


dr] -- — -TQr] 
(—sin^sin 4 >) 


sin( 


dv . dw 


,dv , . dw 

cos^^-sm^^ 


COS(f) 

cos0sin<^|. (C-21) 

We divide the first of these equations by cos^ 6 cos < 56 , the second by cos 6 cos^ (j) and 
the third by cos^ 6 cos^ cf). These divisions assume 6 ^ irf 2 and ^ nf 2. This means 
that the topography cannot have vertical sections along the planes of rotation, i.e. the 
topography function must be single-valued. This is a reasonable assumption given 
that the topography function is assumed to be smooth. After restructuring, 

the three equations become 


tan^ 0 

cos cj) 


+ 1 


du 

dr] 


= tan 0(1 — 


1 


du 


cos (j)) d^ 


— ^tan^ 6 cos + + 


— ^tan^ 0 sin (j) + tan (j>^ 
dv 


dv 


dw tan 6 du tan 9 ^ , 

^ +--tan^sin^— 

d^ cos 9 dr cos 9 or 


tan 9 , ,dw nr • a ^ 1 of a 1 ^ 

-- sm 6 -T —h tan 0 (sm^ — tan<;p) —h tan 0 (cos 9 — Ij —, 

cos 0 ^ drj drj 


(C- 22 ) 


(cos 0 tan^ </) + 1 ) = tan 0 sin 0 -+ (tan 0 + sin 0 tan </>] ^ 

^ dr] cos cp d^ ^ ^ 

.. dw tan 9 du , ( ^ 

an 9 I 

/tan^ (j) 


tdiii (j>du 


dv 




d^ cos <t> dr 
tan 4 ) du 


,cos 0 J dr 


dw 


COS0 

A 

A + 2 /i [ cos (f) ' cos 0 J 
1 / A 1 


-+ tan(?!)(l - COS0) , 

cos 4 > drj dr] 


(C-23) 


tan^ 0 tan^ 6 

+ 


cos (l> \A + 2 /i cos (j) 


11 ^ 

' dr] 
du 


— tan^ $ — — tan 9 tan <f) 


X 


+ 


tan O' 

1 


A 


tan^ <f) 


1 


A + 2/i [ cos 9 cos <f} 


d^ -™-- -r y^^^2/i [cos^ cos^ 
dw tan 9 tan ^ du 
d^ cos 9 cos (p dr 


+ 1 


dv 


1 


+ 

+ 


A 


1 


- tan 4>}^ + 


dv tan <p 


1 


cos 9 y A + 2/i cos 9 
tan 9 / A 1 ^ A 

cos (f> I A + 2 /i cos (j) J dr] 


^ 1 
cos 0 y A + 2/i cos 0 J 

1 dr 


1 tan^ 0 

-) 

A + 2/i 

cos 0 cos (j) 


dv 


(C-24) 


Using the relations tan 61 = dzo{^,T)/d^ and tan<?i = dzo{^,T)/dr cos 9 together with 
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the definitions (38)-(44) leads to the closed set of velocity boundary conditions at a 
free surface topography (35)-(37). 
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FIGURES 


FIG. 1. Rectangular system surface. 

FIG. 2. Curved system surface. 

FIG. 3. The leftmost section shows the topography in a 100 x 100 km area cen¬ 
tered at the NORESS array. The innermost 40 x 40 km was used in the modeling 
experiment. The dotted lines show the positions of the receivers in the two synthetic 
profiles and the circle outlines the NORESS array aperture containing sensors in a 
circular pattern. Labels are in kilometers and elevations are in meters above mean 
sea level. The black area is Lake Mjpsa (123 m above sea level). 

FIG. 4. Snapshots of the vertical component of the particle velocity. The upper 
snapshots are horizontal sections at the free surface and the lower snapshots are ver¬ 
tical sections along a west-east profile through the center of the model. The time 
in the lower left corner of each horizontal section is the time from when the plane 
vertically incident P-wave was reflected from the surface. The straight wavefronts 
parallel to the model boundaries are artificial reflections from the absorbing bound¬ 
aries which are seen also in the vertical sections. Note the strong scattering source 
located at Bronkeberget (2 km N, 10 km E) which is also highly prominent in real 
record analysis (Bannister et ah, 1990). 

FIG. 5. Seismograms extracted at the free surface of the model. The left section 
shows the vertical component of the particle velocity along a 30 km long west-east 
profile through the center of the model. The right section shows 3-component seis¬ 
mograms from a receiver 3.6 km east of the center. 

FIG. 6. Length of topography gradient vector and Frobenius norm of Hessian matrix. 
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FIG. 1. Rectangular system surface. 
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FIG. 3. The leftmost section shows the topography in a 100 X 100 km area centered at 
the NORESS array. The innermost 40 X 40 km was used in the modeling experiment. The 
dotted lines show the positions of the receivers in the two synthetic profiles and the circle 
outlines the NORESS array aperture containing sensors in a circular pattern. Labels are in 
kilometers and elevations are in meters above mean sea level. The black area is Lake Mjpsa 
(123 m above sea level). 


29 







-20 0 20 -20 0 20 

FIG. 4. Snapshots of the vertical component of the particle velocity. The upper snap¬ 
shots are horizontal sections at the free surface and the lower snapshots are vertical sections 


along a west-east profile through the center of the model. The time in the lower left corner 
of each horizontal section is the time from when the plane vertically incident P-wave was 
reflected from the surface. The straight wavefronts parallel to the model boundaries are 
artificial reflections from the absorbing boundaries which are seen also in the vertical sec¬ 
tions. Note the strong scattering source located at Bronkeberget (2 km N, 10 km E) which 
is also highly prominent in real record analysis (Bannister et al., 1990). 
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FIG. 5. Seismograms extracted at the free surface of the model. The left section shows 
the vertical component of the particle velocity along a 30 km long west-east profile through 
the center of the model. The right section shows 3-component seismograms from a receiver 
3.6 km east of the center. 
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FIG. 6. Length of topography gradient vector and Frobenius norm of Hessian matrix. 
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Abstract. We have pursued and compared two- and three-dimensional (3-D) 
finite-difference (F-D) modeling of scattering from free surface topography. A 
velocity-stress formulation of the full elastic wave equations are combined with exact 
boundary conditions for the surface topography and numerically discretized by an 8th 
order F—D method on a staggered grid. We have simulated scattering in 2—D and 3—D 
from teleseismic P-waves using a plane, vertically incident P-wave and real topography 
from a 60 X 60 km area including the NORESS array in south-eastern Norway. Many 
field observations that are not easily explained by simpler 2-D cases are shown to better 
match qualitative effects from 3-D surface topography modeling. These include strong 
amplifications at hills, complex wave pattern caused by scattering, and directivity of 
scattered waves. Snapshots and seismograms show clear conversion from P— to Rg— 
(short period fundamental mode Rayleigh) waves in an area of rough topography in the 
vicinity of the array site. All results are consistent with numerous observations. By 
parallellization of the original software, possibilities have been opened for modeling with 
higher resolution and/or larger areas than before. 
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1. Introduction 

Modeling free surface topography effects is naturally only important for seismic 
profiling on land, whereas in marine settings it satisfies to specify the medium 
parameters without including any explicit boundary conditions. The air-solid Earth 
interface exhibit the strongest possible impedance contrasts. For these reasons alone, 
modeling topography along such a free surface should be important. Any irregularities 
along such a surface would have additional consequences for the results, the more so the 
stronger the gradients and/or irregularities the local topography exhibits [Hestholm and 
Ruud, 1997]. Additionally, uncertain, or even undetermined, parameters are used when 
performing a seismic simulation, like seismic velocities and the associated densities. 
Topography data on the other hand, probably has the smallest error margins of any 
model parameters we use. 

It is a goal for a seismic modeling experiment to achieve good match with 
observations. In this regard the advantage of modeling with surface topography is 
that implicit effects like scattering and conversion will be accounted for automatically 
in the wavefield synthetics. For example, plane incoming P- or S-waves cannot 
convert to Rg-waves when incident on a plane surface of a homogeneous medium. 
Also, P- to SH-phase conversion depends on surface topography. Amplification and 
deamplification of propagating waves can be shown to occur at irregularities and in 
substantial neighborhoods around them [Sanchez-Sesma and Campillo, 1991]. Alluvial 
filled irregular shaped valleys with a plane free surface is also seen to generate strong 
Love and Rg-waves for some incident waves in 3-D [Sanchez-Sesma and Luzon, 1995]. 
Real data correspondence of simulations has proven to be particularly important in 
works on earthquake hazard assessment. In urban areas housing are preferably located 
in hilly, ’socially prestigeous’ areas, often with scant attention to possible amplifications 
from damaging earthquakes. Recently, a study aimed at shedding light on this aspect of 
hazard analysis has been released [Pitarka and Irikura, 1996]. 
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Rg-waves can mask reflections that are the basis of migration. Over the NORESS- 
array in south-eastern Norway the Rg-waves seem to have amplitudes of about 10 % 
of those of the first teleseismic P-wave arrivals [Bannister et aL, 1990; Gupta et aL, 
1993; Hedlin et al.^ 1991]. By quantitatively accounting for topography effects, a dataset 
better suited for migration can be produced. To illustrate how this might be used 
in practice, consider a simulation of an earth model with known surface topography. 

A good correspondence with real data should be obtained to make the initial model 
viable. Then, using the modeling algorithm to propagate the real data backwards in 
time by a technique like inverse time migration, good results might be obtained [Sun 
and McMechanj 1992]. 

With the advent of digital signal recording and deployment of large aperture 
arrays on preCambrian bedrocks in shield areas with moderate topography, unique 
opportunities were given for observing wavefield complexities [Husebye and Ruudj 1989]. 
In case of the NORSAR array (aperture ca 100 km), signal amplitudes were found to 
vary in the extreme by a factor of 10 across the array. In part, this can occur in a 
systematic manner [Haddon and Husebye j 1978]. These observations were modeled in a 
simplistic way in terms of a deepseated (ca 100 km) structural lense producing focussing 
and defocussing across the array aperture. Model residues, apparently of shallow origin, 
were attributed to shallow crust and topography. Another puzzling feature was the 
strong and persistent coda levels in teleseismic recordings. The question was considered 
whether the coda waves were primarily due to source side scattering (S to P and Rg to 
P) or to the receiver side (P to S and P to Rg). This problem was solved by Bannister 
et al. [1990], who identified specific scattering hills in the NORESS siting area where 
the P to Rg conversions were most effective. In array record studies of local events 
P to S scattering appears to be most effective both in forward and backwards modes 
[Charrettej 1991; Hestholm et ahj 1994; Damiy/ 1995]. In profiling surveys scattering " 
effects have been recognized and modeled in 2-D in terms of near surface cavities 
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[Imhof, 1996]. Strong Rg-waves have also been observed and attributed to near-surface 
irregularities with particularly high propagation efficiency in low-velocity weathering 
layers [Levander and Hill, 1985; Ruud et al., 1993]. 

Not much has been published on the modeling of free surface topography in 2-D, 
and even less so in 3—D. A work that considers different slope types and combinations 
of them explicitly in 2-D is that of Jih et al. [1988]. A new 2-D approach [Robertsson, 
1996] incorporates only topography portions parallel to the main axes and classifies 
every surface point in a way similar to Jih et al. [1988]. The scheme can handle both 
elastic and viscoelastic schemes and implements the field variables at each of 7 categories 
of surface topography points separately. A 2—D method which employs a complete 
tensorial formulation of the wave equations for modelling curved interfaces and free 
surface topography is Komatitsch et al. [1996]. This method has the advantage of using 
the same amount of spatial partial derivative calculations as for a Cartesian approach. 
This is 12.5 % less in 2-D and 22.2 % less in 3-D than the chain rule approach used in 
this work. However, it has the distinct disadvantage of 30 % extra memory requirement 
in 2-D and 60 % extra memory requirement in 3-D than the chain rule approach. In 
addition, as in the present method, most chain rule approaches use straight vertical 
grid lines, and therefore the memory difference will be even larger because some partial 
derivatives of the mapping functions vanish. Frankel and Leith’s [1992] 3-D topography 
modeling algorithm employs an F-D scheme of fourth order accuracy in space and uses 
a density taper to zero starting at the level of the free surface. 

Another method for 3-D surface topography modeling of elastic media is the 
boundary integral method in the spatial frequency domain [Bouchon et al., 1996]. The 
medium Green functions have to be calculated for the explicit topography surface as 
integrals over horizontal wavenumbers. The diffracted wave field is the integral over the 
surface topography of the Green functions times unknown source density functions. The 
source density functions are solved for using the conjugate gradient method. Results 
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are shown for an incoming shear-wave polarized along the minor and major axes of a 
cosine-formed elliptically shaped hill. They investigate scattering effects of this hill 
on the wavefield propagating towards an otherwise plane surface in a homogeneous 
medium. Sanchez-Sesma and Campillo [1991] also use a boundary integral method to 
investigate topography effects. Boundary integral methods, or rather their numerical 
discretizations, boundary element methods, have not so far been applied to real 
surface topography modeling. Applications have been restricted to simple geometrical 
structures. This is probably because the discretization applied to each structure must 
be given specific thought. 

The basis of the present method is the vanishing stress condition for a free surface. 
As in Tessmer and Kosloff [1994] (and earlier in Tessmer et al. [1992] in 2-D), a 3-D grid 
is used which is curved in the vertical direction and adapted to the surface topography, 
i. e. the top boundary of the grid coincides with the topography. A coordinate transform 
is used to transfer the elastic, isotropic wave equations from the curved to a rectangular 
grid in which the numerical computations are done. The velocity boundary conditions 
for a free surface are implemented into a local, rotated system at each point of the 
topography. Each of these systems has its vertical coordinate direction coinciding with 
the normal vector of the surface at the given point. The velocity boundary conditions 
are subsequently rotated back to the rectangular system. Once the boundary conditions 
are given in this system, the numerical discretization can be performed. New 3-D 
boundary conditions for the particle velocities at any arbitrary smooth topography 
without vertical subsections have thus been derived [Hestholm and Ruud, 1997]. 

In this study, we state the equations of motion and the surface topography boundary 
conditions for our approach. Next we give a description of the numerical discretization 
and a stability criterium for the method. Then we present simulated scattering from 
teleseismic P-waves using a plane vertically incident P-wave and real topography from 
an area that includes the NORESS array in south-eastern Norway. 
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2. Elastic Wave Modeling Formulation 

The basic equations governing wave propagation in a continuous elastic medium 
are the momentum conservation and the stress-strain relationship. The velocity-stress 
formulation [Achenbach, 1975; Virieux, 1986] can be written in 3-D as in Hestholm and 
Ruud [1997]. We perform a linear transformation from a rectangular (^, r, 7?)-system 
(Figure 1) to a curved (x, y, z)-system (Figure 2), where the relationship between the 
systems is 


x{i,T,ri) = ( 1 ) 

y{^,r,r]) = T, (2) 

2(6d»?) = -^^o(^,r). (3) 

^max 


The topography function zo{^,r) is the local height from the bottom to the surface of 
the curved (x,y, 2 :)-system, and the rectangular {^,T,r])-system is bounded by ^ = 0, 

^ = Uax, r = 0, T = Tmax, T] = 0 and Tj = rjmax- By expanding the velocity-stress 
formulation of the elastic wave equations by the chain rule [Hestholm and Ruud, 1997], 
we get the equations for the wave propagation in the medium 


du 
dv 
dw 
dcF XX 

'W 

dt 

dt 

do'xy 


dcTxx . w>. \do'xx , dcTxy , nft ^ t) —^ + / (4) 






drj 


dt] 


drj 


deXxy A f u \ ^eXxv 


di 

d(Jxz 


+ ^(^,T l) 


'xy 


dr] 

da. 


+ 


dr 

da, 


» + B(f,r.,)% + CK,r)^ + /„ (5) 


drj 

dcTy 


drj 

drxz 


+ + ^ + + C(e,r)^ + /., (6) 


drj 


drj 


drj 


= (A + (|| + ^ (1^ + . (7) 

= A (|| + + (A +2rt (1^ + , (8) 

dw 
drj 


= A(| + 24(e,r,,)| + |^ + S(e,r,,)|j+(A + 2^)C(e,r)^, (9) 




dv 


dv du 


du 




( 10 ) 
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dcTxz 

~w 

dCyz 

dt 


(dw .dw , ^du\ 

(dw .dw xdv\ 

M(^ + BK,T,.l)^ + C’(f,r)^j, 


( 11 ) 

( 12 ) 


where A(^,T,r 7 ), B{^,T,r]) and C'(^,r) are functions of the topography function zo{^,t), 
its spatial derivatives and the rectangular grid coordinates, 


Ci^,r) 


dx zo(^, r) 

^ ^ V dzo{^,r) 
dy zo{^,r) dr 

OtJ ZJmax 

dz Zq(X,t)' 


(13) 

(14) 

(15) 


The topography function 2;o(^, r) is the local height from the bottom to the surface of the 
curved (x, y, 2 r)-system, of which the surface is shown in Figure 2. p is the density and 
A and p are Lame’s parameters, f^, fy and are the components of the body forces, 
and u, V and w are the particle velocity components, a^x, cTyy, <Jzz, <Xxy, (^xz and ayz are 
the stress components. Equations (4)-(12) are the momentum conservation equations 
and Hooke’s law for the curved system, now given in the rectangular (^, r, 77)-system. 

The 3-D particle velocity boundary conditions for a free surface topography given 
in the computational (^, r, 77)-system [Hestholm and Ruud, 1997] can be given by 


<f d'u 
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, dv 


(16) 


(17) 
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with 


and 


+ d ( C 
PiC 


PL 
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^ \ dw ^ I C , 1 

eve"*" 1 dr ^ q \q ) dq 


dw dpdu 1 /( 2 \ 

eq dr e \e ^ ) dr 
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+ pC 
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I 


(18) 


C = 


A -|- Qijji 


(19) 


J dzo{i,T) 

'' = 

e = cos [arctan (d)], 


f = sin [arctan (d)], 
dzoi^,T) 

V = -r- -e. 


q = cos [arctan (p)], 


r = sin [arctan (p)]. 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 


3. Numerical Discretization 

The same spatial and time discretization methods as in the 2-D case are used 
[Hestholm and Ruud, 1994; Kindelan et ah, 1990; Sguazzero et ah, 1989]. The schemes 
employ a staggered discretization stencil of the velocity-stress formulation of the 
elastodynamic wave equations [Levander, 1988] based on Virieux [1986]. Only the 
particle velocities need to be defined at the surface topography. In order to get 
explicit expressions for the particle velocities and stresses at each time step, we stagger 
the vertical velocity component w half a grid length downwards and the horizontal 
velocity components u and v half a grid length positively along the and r-directions 
respectively. Stresses are implemented at the grid nodes or at the middle of the 
’plaquettes’ (the 2-D case is shown in Figure 3). The 3-D boundary conditions 
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(16)-(18) are discretized by 2nd order, staggered F-D operators [Fornberg, 1988], 
from which the particle velocity surface values are solved for from the vertical partial 
derivative expressions. In this procedure, the horizontal partial derivatives are calculated 
one grid length or one and a half grid length below the free surface and considered known 
(Figure 3). Moving from the free surface into the medium, the central, staggered F-D 
method’s order is gradually increased with depth, via 4th and 6th up to 8th order, which 
is the order used inside the domain. The 8th order method is dispersion-bounded and 
cost-optimized [Kindelan et al, 1990]. Along the artificial grid boundaries exponential 
damping is used [Cerjan et al, 1985]. In a layer of 25 points along each grid boundary, 
the stresses and velocities are multiplied by exponentially decreasing terms towards the 
boundary. 

Generally, a stability criterion has to be satisfied in the numerical implementation 
when combining the medium equations (4)-(12) with the boundary conditions (16)-(18). 
It turns out that min{{zo{^,r)} is sufficient for stability in most cases 

[Hestholm and Ruud, 1997]. In such cases, the vertical grid distance (distance between 
grid points) should be a third of the horizontal grid distances in the (^, r, 77)-grid in 
order to maintain the same sampling rate in all directions of the physical (x, y, z)-grid. 

4. Numerical Simulations using Real Topography 

For F-D modeling experiments we have used digital elevation data for an area of 
60 X 60 km containing the NORESS array in south-eastern Norway (Figure 4). This 
hilly area was chosen because of easy access to detailed topography data. Additionally, 
significant P to Rg scattering from specific hills is well documented from NORESS 
record analysis [Bannister et al, 1990]. The code has been parallellized using MPI 
(Message Passing Interface). In this way we are able to run 3-D models with about 
10® gridpoints on the most recent Cray Origin 2000 parallel machine. In the following- 
examples, we display 2-D and 3-D simulations with 0.2 km grid-sampling for a 60 x 
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60 km topography area displayed in Figure 4. The most prominent topography present 
in the dataset is the steep Skreikampen hill immediately west of the southern part 
of lake Mj0sa (the long south-north oriented dark pattern in the south-west area in 
Figure 4). In all the examples shown the incoming wave is a vertically incident plane 
P-wave simulating a teleseismic short period P-phase. The center frequency of the 
Ricker wavelet is 5 Hz, the P-wave velocity of the homogeneous medium is 6.0 km/s, 
the S-wave velocity is 3.46 km/s and the density is 2.0 kgfm^. In this section we show 
simulations for gradually more complex situations, with a 3-D elastic simulation using 
the parallel code as the final result. The outline of displaying results for the same area 
with gradually added complexities will illustrate the need to include these complexities 
in both academic applications and exploration. 

4.1. Topography Surface versus Plane Surface 

Figure 5 shows shapshots for a 2-D simulation of an incident plane wave reflecting 
from a plane, free surface on top of a homogeneous medium of 28.8 km depth. Distances 
are in kilometers, so the area covers a 60 km long horizontal subsection. The left frames 
show the horizontal (li) (top) and vertical (ta) (bottom) particle velocity component 
0.12 seconds after the start of the simulation and 0.04 seconds after the reflection from 
the free surface of the center part of the plane wave. There is 1 second between each of 
the next snapshot frames. The snapshots show a clean reflection from the plane surface 
and artificial reflections from the top corners. All snapshots displayed in this section are 
scaled with respect to a maximum value determined for each simulation. This means 
that amplitudes can be compared directly within the same figure. 

The snapshots in Figure 6 show the exact same times, model depth and source 
location/type as in Figure 5, but here the actual NORESS topography is added (the 
.60 km long, west-east profile of Figure 4). The medium is still homogeneous, so the 
scattering is caused only by the topography. The strongest scattered waves are seen 
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to radiate out from the surface at two points about 25 km and 45 km from the left 
edge, coinciding with areas of steep topography (see Figure 4). The steep hill directly 
west of lake Mj0sa along the profile seems to be at the edge of the present model and 
not contributing much to the scattering because of the absorbing boundary conditions 
acting there. The scattered waves are best seen on the u particle velocity component 
because this component is not affected by the strong vertically oriented wavefront. 
Rg-waves is the dominating feature at the surface, but converted and reflected P- and 
S-waves are also present in the snapshots. The field propagating away from the areas of 
prominent topography has a relatively coherent appearance. 

Figure 7 shows the three following snapshot times for the simulation in Figure 5 
(plane surface). Except for the reflected surface wave seen on the lo-component 
snapshots, only artificial reflections emanate from the top corners and sides. Figure 8 
shows the same snapshots with the topography added, i.e. the simulation of Figure 6 
for the next three times. The artificial reflections propagate further into the medium, 
and the scattered waves from the areas of prominent topography are seen to traverse 
downwards. This is most clearly seen on the w—component. The scattered wavefield 
continues to show a coherent appearance, and no wavefronts are seen after the first two 
clear ones propagating out slightly to the left of the middle of the model. 

4.2. 3-D versus 2—D 

For the 3-D simulation we use the 60 x 60 km NORESS topography area (Figure 4, 
left) for a homogeneous 28.8 km deep model. Exactly the same source position and type 
are used as in the 2-D case, i.e. the same plane, vertically incident Gaussian. 

4.2.1. Snapshots. The snapshots shown in Figure 9 are taken for the same 
particle velocity components (u and w) and at exactly the same locations as the 
topography data for the 2-D case, i.e. an a:^-plane midway along the j/-direction of 
the 60 X 60 km area. The snapshots of the u;-component from the 3-D run seems 
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to differ slightly from the 2-D case by the stronger white band in the center of the 
wavefront. This must be because of an amplitude difference from geometrical spreading 
due to both the numerical discretization and actual physical geometrical spreading. 
Even if analytical geometrical spreading of a plane wave is impossible, it is important 
to remember that after reflection from the surface topography, the wavefront is no 
longer plane. As seen from Figure 9 the scattered surface waves appear to radiate out 
from secondary point sources which coincide with areas of high topographic gradients 
(Figure 4). Energy spreads out into the extra horizontal direction, therefore the 
amplitudes of the u-component, in particular, will be weaker than in 2-D. This is not 
apparent from the figures, since in order to visualize the scattering features, each figure 
is scaled according to a maximum value for that simulation. To the contrary of 2-D 
simulations, in 3-D the strongest amplitudes at a certain time may be from out-of-plane 
scattering, which will be delayed due to its longer travel path and complicate the wave 
pattern. 

Another pronounced feature visible in the 3-D simulation is the disruption of the 
scattered wavefield. It appears to contain much stronger variations in its wavelengths, 
i.e. the spatial spectrum seems to contain more peaks than in the 2-D simulations. The 
total scattered wavefield has a much more complex appearance in 3-D than in 2~D. 
This can be seen on the wavefield pattern for both particle velocity components. 2-D 
simulations are expected to exhibit a greater degree of localization of the scattered 
wavefield than 3-D simulations [Imhof^ 1996]. In view of the apparent energy 
partitioning into frequency peaks in 3-D, this is seen to hold if physical dispersion were 
included (as it would be for viscoelastic wave modeling). 

The increased complexity of the 3-D scattered wavefield pattern is even more 
apparent in the next three snapshots of the u and u;-components shown in Figure 10, 
where the scattered wavefield has propagated into most of-the model. The corresponding 
2-D snapshots (Figure 8) appear quite coherent beside the 3-D ones. Still, the clear 
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wavefronts of the 2-D scattered field can be recognized in 3-D, but they are here more 
complex and interspersed by new wavefronts due to out-of-plane scattering. This 
phenomenon is also the reason for slightly broader wavefronts of the u-component and 
the much more complex general appearance of the 3-D scattered wavefield compared to 
the 2-D one. 

We finally show snapshots for the remaining dimensions from the 3-D simulation. 
Figures 11 and 12 display snapshots of the w particle velocity component for the times 
shown in the previous figures. The upper series in both figures show snapshots along 
the surface topography for the complete 60 x 60 km area, while the lower series is taken 
along the complete 60 x 28.8 km y^-plane at the midpoint along the a;-direction. 

The extra information drawn from the display of the scattered wavefield along 
the surface topography justifies 3-D simulations opposed to 2-D simulations. In 
particular, when investigating surface topography effects, the out-of-plane information 
is invaluable. The positions and shapes of the scattered waves can easily be traced with 
time in the snapshots. The strongest scattered waves can be seen to emanate from 
locations of strong topographic gradients. The strongest effect is seen to be caused 
by the steep area directly west of the Southern part of lake Mj0sa (near the lower 
left area of the horizontal snapshots). This area causes a prominent Rg-wavefront 
that can be traced on all the snapshots. It propagates slowly towards the middle part 
of the topography area, where it can be seen clearly on the last snapshot. Results 
are consistent with the ones found by Pitarka and Irikura [1996] and Bouchon et al. 
[1996], that elevated areas and high gradients amplify wave amplitudes and cause 
scattering away from these areas. The results in Figures 11 and 12 are also consistent 
with 3-D effects found from field experiments and from simulations [Bouchon et aLj 
1996], i.e. the directivity of the wavefield. The most prominent topography causes 
Rg-waves propagating away from it, and these waves have.stronger amplitudes in 
certain directions. The particularly clear wavefront propagating towards the center of 
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the area can be observed in all horizontal snapshots. 

The absorbing boundary condition used in all the simulations is the exponential 
damping method according to Cerjan et al [1985]. Extensive tests have been done 
[Simone and Hestholm^ 1997] to verify that this method is among the best available. 
Nevertheless, a plane incident wave shows the method from its most unfavorable side, 
so that soon after the initiation, reflections will occur and propagate into the domain 
from all grid edges. This is clearly seen from Figures 11 and 12. The reflections appear 
as straight lines of much lower frequencies than the scattered wavefield. They might 
therefore be filtered out from the results by a filtering technique, but one cannot be sure 
of not eliminating some of the physically significant wave portions at the same time. 

The 2-D simulations in this section took about 5 minutes on an IBM RISC 6000 
model 590 workstation. The 3-D run took 2 days and nights (when others were also 
using the machine) on an Intel Paragon machine using 48 processors. This is a much 
slower machine than the 590 workstation, but it was not possible to access enough 
memory on the workstation for this 3-D simulation. However, on a new Cray Origin 
2000 parallel machine using a sequential version of the program on one processor and 
the memory of 8 processors, the time for this 3-D simulation was reduced to 4 hours 
when others were also using the machine. The total memory requirements of the 2-D 
simulation was 2.3 MB, while the size of the 3-D simulation (including all overlap 
domains between processors) was 763 MB. 

4.2.2, Seismograms. This section will compare seismograms for corresponding 
2-D and 3-D simulations. 3-D seismograms are taken from the same 3-D simulation 
that was used for the snapshots. Synthetics are shown along the profile shown as 
the dashed lines in the left part of Figure 4. The first recordings are taken along the 
west-east profile at 40 km from the south grid boundary with stations starting at 25,6 
km and ending at 55 km from the west grid boundary. The second record set is along the 
south-north profile at 40 km from the west grid boundary with the stations distributed 
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from 25,6 km to 55 km from tho south grid boundary. There are 50 stations in each 
profile and the spacing between them is 600 m. In order to compare with 2-D runs, 
we performed 2-D simulations using exactly the same topography data and receiver 
locations as along the profiles for the 3-D simulation. This means that each 2-D profile 
consisted of 50 stations distributed from 25.6 km to 55 km from the lower grid edge 
with an inter-spacing of 600 m. 

Figure 13 compares seismograms for the vertical particle velocity component in 
2-D (left part) and 3-D (right part). The receiver profile is the west-east oriented one 
shown in the left part of Figure 4. The slightly slower 3-D arrival times are due to the 
fact that 4 extra grid points along all edges of the computational domain are included 
in the model compared to the 2-D case. Inside the 3-D domain these 4 points are the 
overlapping layers between processors for the parallel code, but along the computational 
domain edges they are included in the physical model. The strong arrivals striking the 
rightmost receivers and spreading into the model from the right in both seismograms 
are the artificial reflections from the east grid boundary. The other strong arrival 
at all receivers after 5 seconds in the right seismogram of Figure 13 is the artificial 
reflections from the north grid boundary. The left seismogram clearly shows coherent 
wavetrains emanating from the areas of prominent topography. These waves can be 
identified as S- and Rg-waves from their velocities. Going to the 3-D case on the right 
part of Figure 13, the picture becomes much more complicated. Still the prominent 
topographic relief from the 2-D case can be identified as causing the largest amount 
of scattering, but the arrivals are more incoherent and less localized than in 2-D. The 
relative amplitudes of the scattered phases are greatly enhanced compared to the 2—D 
case. As for the snapshots, the 3-D simulation causes the wave pattern to be more 
irregular, which is attributed to out-of-plane topographic scattering. This is also the 
reason for the extra scattered arrivals, in 3-D. The overall weaker-amplitudes of . the 3-D 
simulation (corresponding seismograms are scaled up by 3.5) we attribute to geometrical 
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spreading, which will affect all phases generated near the surface after the plane wave 
has been reflected. In addition, geometrical spreading will occur due to the numerical 
definition of the plane wave at discrete grid points. This effect will therefore weaken all 
3-D arrivals compared to 2-D ones. It is also possible that more destructive interference 
will occur in 3-D than in 2-D because of the extra dimension of discrete point sources 
interacting. 

In Figure 14 we show seismograms from 2-D (left part) and 3-D (right part) 
simulations for the same receivers as in Figure 13, but here the first horizontal particle 
velocity component (the one directed west-east) is displayed. The artificial boundary 
reflections affecting this component, i.e. the ones from the east grid boundary, are seen 
clearly as both P- and S-wavetrains, and after 4 seconds the artificial P-reflections from 
the west boundary can also be seen coming in from the left. The areas of prominent 
topography are seen to give rise to the strongest amplitude amplifications of scattered 
waves from the incident plane wave in both the 2-D and 3-D cases. Still, the appearance 
of the latter is different and more complicated. Additionally, extra scattered arrivals 
can be identified in 3-D. They are caused by out-of-plane topography. Horizontal and 
vertical particle velocity components have the same scaling, and so from comparing the 
two components in each dimension we can confirm that the scattered surface amplitudes 
have the same order of magnitude for all components. 

In Figure 15 seismograms for a 2-D simulation (left) and the 3-D simulation 
(right) are displayed for the 50 south-north oriented receivers shown in the left part 
of Figure 4. The 60 km long 2-D profile with irregular topography is taken along this 
receiver line and covers the complete yz-plane of the 3-D model. The vertical particle 
velocity component is displayed in Figure 15, and 3-D results are again scaled by a 
constant for comparison with the 2-D results. Artificial grid reflections as seen from 
the east boundary for the west-east profile propagate from the north boundary, and 
the large arrival in 3-D after 5 seconds is an artificial reflection from the east grid 
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boundary. In 2-D one can see that the areas of most topographic variations cause 
quite coherent amplitude amplifications. In 3-D the corresponding arrivals get masked 
from scattering from out-of-plane topography, although some of the 2-D topographic 
scattering pattern can be recognized on the 3-D seismogram. The 3-D topography is 
seen to cause scattering of much stronger relative amplitude compared to that caused 
by the 2-D topography. 

Figure 16 shows seismograms for the same 2-D (left) and 3-D (right) simulations 
and for the same receivers as displayed in Figure 15, but here the south-north oriented 
particle velocity component is displayed. The P- and S-reflections from the north and 
south grid boundaries can be identified, as well as strong scattering from the areas of 
prominent topography. Again the 3-D case leads to a much more irregular wave pattern 
and extra scattered waves compared to the 2-D case, although not to the same extent 
as for the vertical particle velocity component. 

5. Discussion 

We have given details on our 3-D F-D scheme for computing synthetic seismograms 
based on the elastodynamic wave equations and produced results pertaining to the 
topography of the NORESS array siting area. The choice of area was motivated by 
easy access to topography data plus many array results exhibiting amplitude variations 
[Haddon and Husebye^ 1978] and scattering phenomena [Bannister et al.^ 1990; Hedlin 
et aLj 1994]. In our synthetic experiments we focussed on two topics; (i) to demonstrate 
the necessity of 3-D modelling for complex topographic features and (ii) scattering 
features tied to prominent hills. We used both snapshots of 2-D displays and profile 
seismograms to illustrate the results obtained in the experiments. Firstly, the 2-D 
synthetics are very much simpler than the corresponding 3-D ones due to out-of-plane 
contributions. This is expected and easily confirmed by f-k or semblance analysis of 
array recordings. A few examples are given by Hestholm et ai [1994]. Here loss of 
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directionality in the Lg-coda waves is demonstrated, i.e. both forward and backward 
scattered wavelets arrive from any direction. Additionally, coda coherency is low or 
hardly significant, at levels of 0.2-0.4 units [Dainty, 1995]. This is in contrast to the 
2-D synthetics snapshots shown in the left parts of Figures 13-16 or from the semblance 
plots in Hestholm et al. [1994]. Also, from Figures 13-16 we see that 3-D coda waves 
are characterized by larger amplitudes than what is the case for 2—D records — naturally 
due to out-of-plane contributions. This reflects a problem in 2-D synthetics, namely 
that excessive and partly unrealistic velocity perturbations (exceeding 5 %) must 
be introduced to match RMS coda levels in observational records [Hestholm et al., 
1994; Dainty, 1995]. 

As observed by Bannister et al. [1990], some hills in the NORESS area literally 
radiate Rg-waves at regular intervals for incoming teleseismic waves of long durations. 
In this regard array recordings are just point observations while the synthetics in 
Figures 6-12 reveal nearly symmetric Rg-radiation from some secondary sources. 
However, with passing time the interference of phases from a multiplicity of secondary 
sources becomes complex (Figure 11), so propagation directionality from a secondary 
source to a receiver will generally weaken. An observational counterpart to this result 
is that Rg-waves rarely propagate further than 60 km in hilly areas typical of the 
NORESS and GERESS (Bavaria, Germany) arrays while across the plains of northern 
Fennoscandia Rg-waves from explosions occasionally propagate out to 600 km. 

We have quoted observations pertaining to strong wavefield amplification on top 
of hills in areas of rough topography. As expected, the synthetics presented here only 
produce moderate amplitude variations along the profiles in Figures 13-16 - less than a 
factor of 2. The reason is that in order to obtain stronger amplitude variations, a large 
slice of the wavefront has to be focussed (compressed) into a small area. This is simply 
not feasible with our homogeneous crustal model. The observations in Haddon and 
Husebye [1978] reflect wavefront ’bending’ taking place more than 100 km away from the 
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free surface sensors. This makes sense with moderate lithospheric velocity perturbations 
on the order of 3-5 %. Undoubtedly, some strong, localized wavefield amplifications 
have taken place in cities destroyed by earthquakes in view of their damage patterns. 

6. Conclusions 

It is demonstrated that 3-D F-D synthetic wavefield analysis including topography 
provides improved insight and a better understanding of surface scattering phenomena. 
Snapshots and seismograms display clear conversion to Rg-waves in areas of rough 
topography. Additionally, 3-D as opposed to 2-D modeling is shown to give significantly 
different results. The added effects from the extra complexities can contribute 
to explanations for deviations between numerical simulations performed on 2-D 
structures and 3-D scattering observed in the field. These discrepancies are connected 
to incoherent and unlocalized scattering, excessive amplitude amplifications and 
directionality dependencies. The last point can be observed in NORESS and GERESS 
recordings. At some sensors the Rg-phase is prominent while hardly visible at nearby 
sensors less than a kilometer away. Many of such wavefield phenomena are out of range 
for 2-D F-D synthetics. The major conclusion from the present study is that there is no 
substitute to 3-D synthetic wavefield modeling for obtaining an adequate comprehension 
of seismic wave propagation in the heterogeneous subsurface. 
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Figure !• Rectangular system surface. 



Figure 2. Curved system surface. 



Figure 3. Computational staggered grid, (^, T 7 )-plane. 
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Figure 4. Left: Map shows the topography of the 60 x 60 km area used in the 3-D 
simulation. The dashed lines show the positions of the receivers in two profiles and the 
circle outlines the NORESS array. Labels are in kilometers and elevations are in meters 
above mean sea level. The black area to the south-west is Lake Mj0sa (123 m above sea 
level). Right: Topography profiles along the two lines each of 60 km length shown on the 
map. They are midway along the y-direction and x-direction and cover respectively the 
complete x- and y-dimension of the area. Horizontal axes are in kilometers and vertical 
axes are in meters above mean sea level. 
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Figure 5. Snapshots of 2 D simulation of a plane, vertically incident P—wave reflecting 
from a plane, free surface. The first snapshot shows the wavefield 0.04 seconds after 
its reflection from the surface, and the time lap between each of the next snapshots is 1 
second. Distances are in kilometers, and the upper and lower series display the horizontal 
and vertical particle velocity components respectively. 
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Figure 6. Same as Figure 5, but for an irregular free surface boundary condition. The 
topography is that of the west-east profile of Figure 4. 
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Figure 7. Same as Figure 5 for the next three snapshots. 
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Figure 8. Same as Figure 6 for the next three snapshots. 
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Figure 9. Snapshots of 3-D simulation of a plane, vertically incident Ricker wavelet 
reflecting from an irregular free surface. The first snapshot shows the wavefield 0.04 
seconds after its reflection from the surface, and the time lap between each of the next 
snapshots is 1 second. Distances are in kilometers, and the upper and lower series display 
a horizontal (u) and vertical particle velocity component respectively. Snapshots are 
shown at the midpoint of the direction of the second horizontal dimension and correspond 
exactly to the snapshot locations of Figure 6 for the 2-D case. 
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Figure 10. Same as Figure 9 for the next three snapshots. 
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Figure 11. The simulation displayed in Figure 9 for the same times, but here snapshots 
are shown at the surface topography and along the plane of the midpoint of the first 
horizontal direction (along the yz-plane). Both series display the vertical particle velocity 
component (in). 
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Figure 12. Same as Figure 11 for the next three snapshots. 
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Figure 13. Left: Seismogram from 2-D simulation of a 60 km long topography data 
profile taken along the west-east oriented receiver line from the left part of Figure 4. 
The profile is 28.8 km deep and extends from 25.6 km west of the leftmost receiver to 5 
km east of the rightmost receiver. The profile contains the displayed line of 50 west-east 
oriented receivers inter-spaced by 600 m. The vertical particle velocity component is 
displayed. Right: Seismogram from the 3-D simulation of the 60 x 60 x 28.8 km model 
displayed in the snapshots. The 50 receivers are located along the west—east oriented 
receiver line displayed on the left part of Figure 4 and inter—spaced by 600 m. Because 
of geometrical spreading, the amplitudes are multiplied by 3.5 to make the main arrivals 
the same order of magnitude as those for the 2-D case on the left. The vertical particle 


velocity component is displayed. 


























































Figure 15. Left: Seismogram from 2-D simulation of a 60 km long topography data 
profile taken along the south-north oriented receiver line from the left part of Figure 4. 
The profile is 28.8 km deep and extends from 25.6 km south of the receiver furthest to 
the south to 5 km north of the receiver furthest to the north. The profile contains the 
displayed line of 50 south-north oriented receivers inter-spaced by 600 m. The vertical 
particle velocity component is displayed. Right: Seismogram from the 3-D simulation of 
the 60 X 60 X 28.8 km model displayed in the snapshots. The 50 receivers are located 
along the south—north oriented receiver line displayed on the left part of Figure 4 and 
inter-spaced by 600 m. Because of geometrical spreading, the amplitudes axe multiplied 
by 3.5 to make the main arrivals the same order of magnitude as those for the 2-D case 


on the left. The vertical particle velocity component is displayed. 
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Figure 16. Seismograms from the same 50 receivers and the same 2-D (left) and 3-D 


(right) simulations as shown in Figure 15, but here analogous horizontal particle velocity 
components (the south-north oriented ones) are displayed. 














































Appendix 8 




1 




RECOGNIZING EXPLOSION SITES WITHOUT SEISMOGRAM READINGS: 
NEURAL NETWORK ANALYSIS OF ENVELOPE TRANSFORMED 
MULTISTATION SP RECORDINGS 3 - 6 HZ. 


Yu. V. Fedorenko^’^^ E. S. Husebye^^ B. Heincke^’^^ and B. O. Ruud^^ 

1) Institute Solid Earth Physics, University of Bergen, Norway 

2) Institute of North Ecology Problem, Fersman street 14, Apatity, Russia 

3) Institute of Geophysics, Albrecht University, Kiel, Germany 

Received 1997 November 28; in original form 1997 August 15. 

SUMMARY 

Seismic waveform similarities for closely spaced earthquakes and explosions in 
particular are well established observationally. In many industrialized countries of low 
seismicity more than 90 % of seismic event recordings stem from chemical explosions and thus 
contribute significantly to the daily analyst workload. In this study we explore the possibility of 
using envelope waveforms from a priori known explosion sites (learning) for recognizing 
subsequent explosions from the same site excluding any analyst interference. To ensure high 
signal correlation while retaining good SNRs we used envelope transformed waveforms 
including both the P and Lg arrivals. To ensure good spatial resolution we used multistation 
(network) recordings. The interpolation and approximation neural network (lANN) of Winston 
(1993) was used for learning the computer to recognize new explosion recordings from a 
specific site using detector output event files of waveforms only. The lANN output is a single 
number between 0-1 and on this scale an acceptance threshold of 0.4 proved appropriate. We 
obtained 100 % correct decisions between two sets of 'site explosions' and hundreds of ‘non¬ 
site’ explosions/earthquakes using data files from the Norwegian Seismograph Network. 
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INTRODUCTION 

An experienced analyst is often capable to locate an earthquake simply by a brief visual 
inspection of a single station recording (Vesanen, 1944). With easy access to digital recording 
the ability to recognize specific explosion sites at least in a semi-automatic manner has been 
explored by several investigators like Israelsson (1990), Harris (1991), Riviere-Barbier and 
Grant (1993) and lost, Schweitzer and Harjes (1996) among others. The basic idea is that 
many explosions are confined to a few sites which should be recognizable through correlation 
and cluster analysis of waveform records. Research focus has been mainly on separating closely 
spaced explosions, although accurate “ground truth” information is often lacking. Since the 
source correlation distance may exceed a few wavelengths of the dominant P-signal period it is 
not easy to resolve closely spaced explosion sites. To improve resolution, ‘covariance records 
in combination with clustering analysis and a manual touch was used by Israelsson (1990) and 
Riviere-Barbier and Grant (1993), while Harris (1991) expressed site resolution in a 
probabilistic manner in terms of bandwidth and time duration. Joswig (1990) on the other hand 
used sonograms for seismic signal detection based on a priori known signal patterns. However, 
the above schemes do not appear convenient for semi-automatic seismogram analysis of array 
and network recordings - not easy to use nor robust enough. Nevertheless, the logic of using 
seismic waveforms of various kinds for identifying local explosion sites is attractive since it 
would ensure to pinpoint locations with a minimum of analyst effort. Such a scheme would be 
of particular interest in low-seismicity countries where up to 95 % of local event recordings are 
explosions. In this short note we describe a neural network scheme for recognizing specific 
explosion sites using envelope transformed (short period vertical component) network 
recordings of local events. 
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DATA ANALYSIS SCHEME 

The problem considered here is a variant of classical seismic source discrimination ~ 
earthquake or explosion (Husebye and Dainty, 1996). However, it is simpler, as we aim to 
recognize explosion sites using multistation envelope waveform patterns tied to the basic 
assumption that these are site specific. These waveforms are envelope transforms of bandpass 
filtered (3-6 Hz) recordings, formed by first computing signal power in 2 sec long moving 
timewindow (STA window) at 0.5 sec intervals. These STA time series are futher smoothed by 
convolution with 3 ‘boxcar’ functions of length 6, 8 and 10 sec. Finally, the square root of 
each sample in the new time series are computed so the envelopes effectively consist of 
smoothed RMS amplitudes. 

Various kinds of neural network schemes have proved to be efficient tools in 
discrimination analysis (Dowla, 1996). For the site recognition problem the so-called 
interpolation and approximation neural nets (lANN; Winston, 1993) or radial basis function 
networks (Dowla, 1996) appear to be ideally suited for our purpose. It has the virtue that 
training (or learning) can be tied to solving a set of linear equations hence in this case it is not 
necessary to involve any optimization method in order to obtain the weights. Below we 
outline the construction of this particular variant of neural nets following the approach of 
Winston. 

The starting point is several inputs x = (xj, X 2 ,...,x^) und one output Y as 

illustrated in Fig. 1. The input vector x consists of K sub-vectors, x —(xj, X 2 ,...,x^) 
where K is the total number of stations used in analysis. Each sub-vector x„, m = l,2,K 
is the envelop record for station m. The number of samples in each sub-vector is different 
because it is determined using the distance from each station to the site when calculating time 



4 


windows. Lg dominates our waveforms, thus the different station traces are weighted 
accordingly. Station weighting, whether tied to body or surface waves geometrical spreading 
was not critical for the final results. The output T may be a single number (the score) and 

preferably Y e [o, l] for facilitating performance comparisons. The aim is to predict future Y- 
values for new (and unknown epicenters) given a proper data set (or learning set) of input- 
output combinations. 

Let the output function T(x) have the following properties; (i) the value of Y is 
exactly equal to 1 whenever the input event coincides to the one of the events in the learning 
set. (ii) The value of Y should be close to 1 for all other events stemming from the site subject 
to recognition analysis and close to 0 for the remaining events. The output Y is denoted the 
'interpolation' function by Winston (1993) and may be expanded into a weighted sum of other 

functions /;.(x), Y(x) = Y,l^WJ,ix), where 5 is the number of events in learning set. It 

consists of multistation records from the specific site we want to recognize. 

Using vector notation this might also be written as 

F(x) = W^f ( 1 ) 

where W = (w,, ,.., f is a vector of weights. 

According to Winston, a convenient choice for /((x) is a Gaussian function that may 

be written in the form /i (x) = exp|-'^||x —c, || j, where the norm ||x —C;|| 

and N is the number of observation points (total number of samples) in input. The shape of 

y.(x) =/(|||x-c,l|^,crj is controlled by given values of <7 and C;. The latter is a reference 

vector obtained from the i -th event of the learning set. Because each /;(x) -value depends on 
only one c, -value the /;.(x) -function together with the weight are tailoring the influence of 
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the i -th learning sample on future predictions of F. To summarize, the function Fis computed 
by a 2-layer net by which the single node in the second layer (see Fig. 1) computes a weighted 
sum of the outputs of the first layer nodes. In turn, each of the first layer nodes “computes a 

Gaussian function centered on the “sample” input. 

For a given o - value of /, a set of linear equations can be written for calculating the 

vector W in eq. (1). In our lANN approach the output F is 1 for samples drawn from the 
learning set, so, substituting the input x by the c,. we rewrite equation (1) 


X/yc,-c,|| ,cTjlI^ = l,for 7 = l,..,S (2) 

or, in matrix notation AW = 1 where l = (l,l,...,l) , Ay =/( 


Xhe linear system of eq. in (2) comprises S unknowns but also S- 

equations. The A-matrix is positive definite and symmetric. However it may be ill posed due to 
the arbitrary choice of events for the learning set (some very small eigenvalues may appear if 
chosen envelops are quite similar) so the SVD approach (singular value decomposition) was 
preferred for solving eq. (2). The record lengths vary between individual stations. Preferably 
both the P and Lg envelopes should be included for all stations in a learning set, which in turn 
should comprise the largest number of network recording stations available because this 
“master” learning set then can easily be adjusted to any sub configuration of the network. This 
property is important since quite often data from one or more stations in a network are missing 
and corresponding new weights are fast to compute from eq (2). Using network envelopes has 
three advantages, (0 high signal similarity, so few events are needed in learning set; 
(ii) robustness, since envelope records are used; and (iii) quite good spatial resolution; since we 
use recordings including both P- and S (Lg)-waves frorn widely spaced network stations. Note, 
using high-frequency (8-16 Hz) explosion records from a single array, Israelsson (1990) and 
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Harris (1990) obtained 2-4 km spatial resolution while ours is around 6-8 km (Fig. 5). 
However, as demonstrated below, envelope processing is far more robust and hence 
convenient for automated record analysis. 

DATA ANALYSIS AND RESULTS 

The data used are local event recordings from the Norwegian Seismograph Network 
(NSN) (Fig. 2) for the period September 1992 to May 1997. The network is operated in a 
trigger mode implying that original waveforms from individual detecting stations are merged 
into easily accessible event files. An essential element in our site recognition analysis is to 
estabUsh the reliable learning set. In this regard, only 2 explosion sites proved convenient, 
namely the open pit mine Titania (TITA) and an underwater harbor construction site in 
Geiranger Fjord (GEIR). For both of these sites written notifications of each shot were 
received. The other events subjected to analysis were different for TITA and GEIR reflecting 
varying station detectabilities for these 2 sites (Fig. 2). Since P- and S-group velocities are 
quite different the envelope shapes change rapidly with epicenter distance thus ensuring 
adequate spatial resolution. Origin time is not a problem since for the a priori known sites we 
can accurately calculate/predict the relative move-out times between the station used. 
Alternatively, one may use one event in the learning set to set the time windows manually for 
each station without any knowledge of the velocity model or station distance. Through ‘sliding 
window operation’ the lANN recognition scheme amounts to a two-pronged test; (i) that 
signal record move-out times correspond to the specific sites and (ii) their envelope waveform 
shapes conform to those of the learning set (Fig. 3). 

The envelope event files were screened for appropriate candidates conditioned on; (i) 
relatively large number of recording stations (3-7) and (ii) a minimum of 3 stations in common 
with the learning set. On average, the Titania and Geiranger explosions are recorded at 3 - 6 
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stations only so there is an obvious need to adjust the learning set to the specific event 
recordings at hand as discussed above. Even in the extreme case of only 2 recording stations, 
site recognition was sometimes feasible. In Fig. 3 the Titania and Geiranger learning set 
envelopes are shown. The site recognition performance depends on the variance (<T- 
parameter) in the presumed Gaussian distribution (1, 2). Testing gave that the best population 
separation was obtained for a cr value at 0.04 which was used in subsequent analysis. 

The outcome of the site recognition analysis are displayed in Fig. 4 for Geiranger and 
for both sites in Fig. 5. The performance for Geiranger was better than that for Titania which 
partly is attributed to a better network configuration (Fig. 2). Occasionally, “true” Titania 
events had scores below the threshold of 0.4. To our surprise this was due to small errors of a 
few seconds in the station timing system (ODDI, ASK, EGD). The site recognition scheme 
also has good spatial resolution. As apparent fi-om (Fig. 5), any event more than 6-8 km away 
from ‘true’ site location would be rejected. If closer and with similar envelope shape, it could 
possibly be accepted as a ‘site explosion’. Best results were obtained for Geiranger where with 
one exception aU scores for confirmed explosions exceeded 0.7 units. The non-Geiranger event 
scores were mostly zero and never exceeded 0.2 units. For Titama the results are sunilar but 
somewhat lower in absolute terms, i.e., score between 0.4 and 0.7 are relativly frequent. These 
lower scores are in part attributed to the aforementioned timing errors. However, some non- 
Titania events in the range 100 - 150 km from Titania have scores between 0.2 - 0.4 and the 
corresponding epicenter locations are in some cases close to the stations KMY or BLS t k 
(Fig. 2a). The explanation here is as follows: only 3 stations are used namely ASK, EGD and 
KMY or BLS with an appearant lined-up versus the Titania .site. .In two cases, low site 
recognition scores were due to interference from other events causing Lg envelope shape 


modifications. 
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DISCUSSION. 

Our novel approach to explosion site recognition has, unsurprisingly, performed well 
on this dataset. The discrimination power of our lANN approach is surprising; rather small 
timin g errors either stemming from clock defects or other epicenters a few kilometers away 
from site suffice for rejecting the event in question. This also implies that our approach as 
presently formulated hardly would be applicable for earthquake record analysis because of less 
spatial clustering of foci even in cases of swarm activities. 

In areas like Fennoscandia and large parts of Europe, site recognition scheme has the 
potential of significantly reducing tedious analyst effort in bulletin preparations since 90 - 95 % 
of all local events stem from mining and construction explosion activities. The former seismic 
source type is most convenient to handle since operational span often amounts to many years 
while construction works finish in months or at most a few years. The latter can also be 
handled effectively by our selection approach since the learning set can be established with 5 - 
10 events. However, a potential drawback is that analyst is often unaware of temporary 
explosion sites, leading to ground truth problems in event selection for learning sets. Also 
bulletin epicenter solutions are quite error prone (Fig. 2a and 2b) so this parameter is not 
helpful in this regard (Husebye, Ruud and Dainty, 1998). 

Our site recognition scheme does not incorporate specific earthquake/explosion source 
type criteria but essentially discriminates on spatial locations for presumed identical envelope 
recordings. However, for this type of problems we may for envelope generation use several 
bandpass filters reflecting those Pn and Lg signal frequencies for which the best 
earthquake/explosion source separation has been obtained (Dowla, 1996). In special cases, we 
may consider non-envelope transformed record analysis like those of Israelsson, (1990), 
Joswig (1990) and Harris, (1990). Regrettably, to obtain an adequate data base for recognizing 
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source type appears to be problematic. The main potential of our flexible explosion site 
recognition approach is to reduce significantly the analyst work in many industrial countries. 
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FIGURE CAPTIONS 

Figure 1. The Interpolation and Approximation Neural Net (lANN, Winston, 1993). 

Figure 2a. The stations of the Norwegian Seismograph Network (NSN) used in our site 
recognition analysis for the open pit mine explosions at Titania (TITA). Rings are event 
epicenter locations taken from the NSN-buUetins while ring shading indicates recogmtion 
score: white = 0.0-0.1; gray = 0.11 - 0.39 and black > 0.40. The dispersion of high score 
events around the Titania mining site simply reflects errors in the local bulletin for these 
explosions. Such errors are significantly smaller for the Geiranger explosions. Score results 
also presented in Figs. 4 and 5. 

Figure 2b. The stations of the Norwegian Seismograph Network (NSN) used in our site 
recognition analysis for the underwater construction explosions in Geiranger. In comparison to 
Titania (2a) bulletin location errors are small. Figure caption otherwise as for Fig. 2a. 

Figure 3a. Stippled envelope traces for individual station recordings for known (confirmed) 
Titania explosions which were used as learning set in the site recognition analysis. These 
recordings, typically comprising 3-10 events, are flexible in the sense that number of stations 
used can be reduced to 3 or even 2 if some of the station recordings are lacking or defect. 
Notice the consistent dominance of Lg-phase and the often weak, first arriving, Pg-phase. The 
latter is occasionally lost in the ambient noise. For convenience, traces are time shifted relative 
to expected P-onset time computed from epicenter distances. The traces marked by a solid line 
are a ‘new’ event subject to lANN recognition test. The score is given on the top of the figure 
together with the c -value of 0.04 and Origin Time (year, day, hr, min, sec). 

Figure 3b. Stippled envelope traces for individual station recordings for known (confirmed) 
Geiranger explosions which were used as learning sets in the site recognition analysis. Caption 
otherwise as for Fig. 3 a. 
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Figure 4. Geiranger site recognition scores (vertical axis) as a function of distance between the 
site and the epicenters depicted in Fig. 2b. Regrettably, we lack information about which 
events are explosions and which are earthquakes. Notice that the scores, with one exception, 
are uniformly high, i.e., larger than 0.7 units. The non-site (more than 600) events have 
unif orm low scores below 0.2 units so the discrimination between site and non-site events is 
excellent. 

Figure 5a. Simulated Titania score performances as a function of fictitious epicenter position 
relative to the 'true' site location. The recognition score acceptance of 0.4 units is equivalent to 
an area of radius 8 km. 

Figure 5b. Simulated Geiranger score performances as a function of fictitious epicenter 
position relative to the 'true' site location. Due to a non-uniform network configuration the 
recognition score acceptance area is elliptic with the largest half-axis of 10 km. 
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4.9 Signal site recognition using single station 3-component records 

As mentioned, seismic waveform similarities for closely spaced explosions are well established ob- 
servationally. Also, encouraged by our previous site recognition results we considered additional 
flexibility in our approach to this problem essentially how to incorporate additional seismic record 
information in the neural network analysis. Another motivation here was to safeguard against 
technical defects like timing errors which otherwise would ruin multistation analysis. The alter¬ 
native strategy choosen for extracting relative comprehensive signal attributes was that of using 
single station 3-component (3C) records. By introducing complex demodulated record operations 
(including both the P- and S-wave arrivals and part of the coda) the 3C covariance matrix would 
contain 3x3 elements of different time information. From other source classification studies we 
know that both explosion and earthquake spesific information are contained in several frequency 
bands so we formed covariance matrixes for a suit of 12 different frequency bands. In other words, 
a single 3C station event record produced 9x12 = 108 pieces of time information- Having access 
to ’ground truth’ information from 2 underwater construction sites (Mongstad and Geiranger) 
we had 4-6 events (enough) for site recognition learning sets. Obtained results were excellent; 
all events from the construction sites (not in learning sets) were classified as such while all other 
events (more than 200) were classified cis non-site events whether source was an explosion or an 
earthquake. Presuming log-normal distribution of the neural network output score factor we esti¬ 
mated that the probabilities of false alarm or missed ’detection’ were less than 1:1,000,000. Also, 
spatial resolution was good (less than 10 km) even perpendicularly to the station-site azimuth 
direction. During this 3C site recognition exercise we found in the epicenter map used, 2 cluster of 
events with which origin times were mainly between 10, - 14. hours that is prime working hours. 
Selecting events directly from bulletin for learning sets but without ground truth information at 
hand did produce results somewhat less convincing in comparision to those obtained for Mongstad 
and Geiranger. This kind of problems, lack of ground truth information for appearant event clus¬ 
ters, would in future be addressed through an initial multivariate analysis scheme in order to sort 
out ’most similar’ event recordings and hence the likely event candidates for learning sets. 
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